# Nonparametric monitoring of sunspot number observations: a case study

Sophie Mathieu
ISBA/LIDAM, UCLouvain
and
Laure Lefèvre
Solar Physics and Space Weather department, Royal Observatory of Belgium
and
Rainer von Sachs
ISBA/LIDAM, UCLouvain
and
Véronique Delouille
Solar Physics and Space Weather department, Royal Observatory of Belgium
and
Christian Ritter
ISBA/LIDAM, UCLouvain
and
Frédéric Clette
Solar Physics and Space Weather department, Royal Observatory of Belgium

###### Abstract

Solar activity is an important driver of long-term climate trends and must be
accounted for in climate models. Unfortunately, direct measurements of this
quantity over long periods do not exist. The only observation related to solar
activity whose records reach back to the seventeenth century are sunspots.
Surprisingly, determining the number of sunspots consistently over time has
remained until today a challenging statistical problem. It arises from the
need of consolidating data from multiple observing stations around the world
in a context of low signal-to-noise ratios, non-stationarity, missing data,
non-standard distributions and many kinds of errors. The data from some
stations experience therefore severe and various deviations over time.
In this paper, we propose the first systematic and thorough statistical
approach for monitoring these complex and important series. It consists of
three steps essential for successful treatment of the data: smoothing on
multiple timescales, monitoring using block bootstrap calibrated CUSUM charts
and classifying of out-of-control situations by support vector techniques.
This approach allows us to detect a wide range of anomalies (such as sudden
jumps or more progressive drifts), unseen in previous analyses. It helps us to
identify the causes of major deviations, which are often observer or equipment
related. Their detection and identification will contribute to improve future
observations. Their elimination or correction in past data will lead to a more
precise reconstruction of the world reference index for solar activity: the
International Sunspot Number.

Keywords: Statistical process control; Support vector machine; Correlation;
Missing data; Control chart; Block bootstrap

##  1 Introduction

The main objective of this paper is to provide an efficient monitoring for the
sunspot numbers. Until recently, such a monitoring could not be developed with
the available statistical methods due to the complexity of the series.
Moreover, existing modeling and error quantification of the data were lacking.
Advances in statistical process control for panel data combined with a
thorough uncertainty model for the sunspot numbers (Mathieu et al.,, 2019)
allow us to construct an appropriate monitoring method for these data.
Our developed method bears however the potential to be applied to a much wider
range of problems. Hence, this paper follows two threads, one related to the
specific patterns and considerations of sunspot observations, the other to the
development, selection, and parametrization of a statistical methodology for
monitoring them.

###  1.1 Overview

Sunspots appear as dark areas on the Sun. They correspond to regions with a
high local magnetic field and relate to the solar activity. The International
Sunspot Number is an indicator of this activity. It is composed of the number
of sunspots, $N_{s}$ and the number of sunspot groups, $N_{g}$ via a
composite, $N_{c}=N_{s}+10N_{g}$. This index is one of the most intensely used
time-series in astrophysics (Hathaway,, 2010). It enters into models of the
influence of the solar irradiance on the Earth climate (Haigh,, 2002; Ermolli
et al.,, 2013) and in space weather predictions (Temmer et al.,, 2001; Wang
and Colaninno,, 2014).

Although astronomers started observing sunspots in the beginning of the
seventeenth century, it remains surprisingly difficult to arrive at an
accurate daily determination of their numbers. Three difficulties stand out:
observability, resolution and interpretation. For Earth-based observatories,
the Sun cannot be observed when there are clouds. Instruments with different
resolutions may give rise to different counts of sunspots. Distinguishing
sunspots and groups of sunspots require experience and even experts sometimes
disagree. The estimate is therefore a weighted consensus among all
observatories, also called “stations” participating in the effort. Different
observers may vary in skill and their skills may vary over time. Some stations
observe systematically more sunspots than others and there may be shifts due
to instrument changes. Moreover, the observed series from different stations
cover different time periods. We are therefore facing panel data with missing
values of various patterns, time-varying biases and non-overlapping time
periods.
There are also intrinsic sources of variability. Some sunspots are only
visible over short periods. Their number may therefore change during the day.
Finally, the sunspot activity itself is subject to substantial variability at
multiple time-scales.
The statistical challenge is to untangle observation-related variability from
actual variations of the solar activity: to detect and correct the former and
to extract the latter. These variability patterns are non-normal and
autocorrelation naturally occurs.

There are also many disturbances such as jumps, drifts or oscillations. Those
have been partially studied on the short-term in the previous works by Morfill
et al., (1991), Vigouroux and Delache, (1994) and later Dudok de Wit et al.,
(2016). More recently, Mathieu et al., (2019) developed a comprehensive
uncertainty model that involves three types of station-dependent errors, at
short-term, long-term time scales as well as during solar minima. It reveals
the various irregularities that stations may exhibit.
Assuring the quality of the consolidated series therefore calls for an
automated tool for supervising the observations in quasi real-time. This
procedure should monitor the stations and send alerts when they start
deviating to prevent the occurrence of large drifts. Owing to methodological
advances in the sunspot numbers uncertainty modeling and in statistical
process control (SPC), it is now possible to develop such a method.

###  1.2 Related work

Different procedures have been developed in the SPC literature to monitor
complex panel of data. The aim of this case-study is not to compare them. It
is to adapt the most-promising existing procedure to efficiently control the
deviations of the sunspot numbers over time. Therefore, we propose in the
following a method based the dynamic screening system developed by Qiu and
Xiang, (2014).
The method of Qiu and Xiang, (2014) is based on extensions of the classical
CUSUM (Page,, 1961) chart. It is composed of two steps. First, the regular
patterns (i.e. the mean and the variance) of the data are estimated on a
subset of non-deviating or in-control (IC) series. Second, the data are
standardized by these patterns and monitored by a CUSUM chart designed by a
block bootstrap method. This procedure constructs, without any parametric
assumption, a control scheme that is valid for non-normally distributed and
serially correlated data. Moreover, the regular patterns of the data are
estimated over a subset of stable series. This allows the chart to detect
shifts in the mean level of each series, where the series is allowed to have a
mean changing over time. Hence, the method can accommodate the quasi-periodic
variations in the sunspot number series which are related to the eleven year
solar cycle (Hathaway,, 2010). This cycle is intrinsic to the signal and
unlinked to the long-term stability of the observations. The control scheme of
Qiu and Xiang, (2014) seems therefore particularly appropriate for our
problem. It allows us to treat the sunspot data which are non-stationary, non-
normally distributed and autocorrelated.

In the following, we use and extend the work of Qiu and Xiang, (2014), to
bridge the gaps between the method and the specific requirements of our
problem. Those gaps are two-fold. (1) The method of Qiu and Xiang, (2014) —as
all other methods that we encountered in the literature— cannot be used
however without knowing a priori which stations are in-control. This
information is not available for our data, where _all_ stations are expected
to contain several kinds of deviations (jumps, oscillating shifts, etc) in
their observation period. (2) The method operates with a control chart which
sends an alert when a deviation is detected, yet without providing any
information about the nature of the shift. Such information is however crucial
for us, since it allows to further investigate the causes of the shifts.
Although several methods have been developed to automatically predict the size
of a shift after an alert (see for instance Cheng et al., (2011) and the
references therein), they are not adapted to data which are (simultaneously)
non-normally distributed, serially correlated and contaminated by strong
noise.

###  1.3 Aims

In the following, we propose a nonparametric monitoring that is tailored to
the complex features of the sunspot numbers: (a) the missing values, (b) the
strong noise, (c) the complex autocorrelation structure and (d) the non-
normality. Our method extensively exploits the information contained in the
panel to establish a robust IC reference from the network. This allows us to
monitor the stations without prior information on their stability. We complete
the method by a support vector machine (SVM) procedure that efficiently
predicts the size and the shape of a shift once an alert has been raised.
Although we could manually build a library with typical shapes and sizes to be
compared to the deviations, we select the automatic SVM approach instead.
The control scheme is then applied on past observations to study the
deviations of the sunspot numbers. The procedure automatically detects major
deviations identified recently by hand in some stations. It also unravels many
other deviations, unseen in previous analyzes. In particular, small and
persistent shifts that are difficult to identify manually are detected by the
method. The precise information about the deviations predicted by the SVM
procedures allows us to determine the causes of some prominent deviations.
This sets the ground for a future enhancement of the quality of the series.
Moreover, the monitoring procedure provides the possibility to be used in
real-time to preserve the long-term stability of the stations. It also paves
the way to a future redefinition of the International Sunspot Number based on
several stations that are stable over time.

This article is structured as follows. In Section 2, we present the main
properties of the data and their model. The methods are explained in Section
3. This includes the complete monitoring scheme as well as the SVM procedures
to predict the size and shape of the deviations. In Section 4, we apply the
proposed method on the sunspot data at different scales, in order to detect
both high- and low- frequency shifts and discuss the results on actual
stations. In a final section 5, we give some concluding remarks and
perspectives. Supplementary materials provide some more details about the
monitoring scheme as well as more examples of monitored stations.

##  2 Data

The data and their specific features are first presented in this section.
Then, we introduce the uncertainty model associated to the sunspot numbers.
The component of the model that will be monitored in this paper is finally
presented alongside with its estimating procedure.

###  2.1 Presentation of the dataset

The period under study embraces the most recent part of the series and extends
from January 1, 1981 (when the sunspot numbers production center moved from
Zurich to Brussels) till December 31, 2019. It ranges from the descending
phase of solar cycle (SC) 21 to the minimum between SC 24 and
25111https://en.wikipedia.org/wiki/List_of_solar_cycles and covers thus three
complete solar cycles. The data are composed of the daily observations of a
network of 278 Earth-ground observatories disseminated across the world. The
records contain the number of spots $N_{s}$, groups $N_{g}$ and composite
$N_{c}$. They are distributed through the World Data Center Sunspots Index and
Long-term Solar Observations (WDC-SILSO)222 The data are available at the
following link: http://www.sidc.be/silso/. In the following, we denote by $t$,
$t\in 1,...,T$, the date-time of the observation and represent the index of
the stations by $i$, $i\in 1,...,N=278$.

The data have complex features that are described below. Those should be taken
into account in the design of the monitoring procedure.

  1. 1.

As studied in Mathieu et al., (2019) and previous works, the data are by
nature non-normally distributed.

  2. 2.

The data contain around 70% of missing values over the period studied. Those
are mainly caused by the non-overlapping observing periods of the stations.
Indeed, some stations started observing only recently while older stations
stopped their activity well before 2019. The stations also contain various
percentages of missing values (ranging from 15% to 75%) over their active
observing period, mainly due to weather conditions that prevent the
observation of the Sun.

  3. 3.

The stations are also correlated across the panel and along time since a
sunspot may stay from several minutes up to several months on the solar
surface.

  4. 4.

Due to the observing conditions and the solar variability, the series
experience a wide range of deviations that vary in shape and size.

###  2.2 Model

Let $Y_{i}(t)$ represent either the number of spots, groups or composite
observed in station $i$ at time $t$. The observations may be decomposed into a
common solar signal, generically denoted by $s(t)$, corrupted by three types
of station-dependent errors in a _multiplicative_ framework (Mathieu et al.,,
2019):

$\begin{split}Y_{i}(t)=\left\\{\begin{array}[]{ll}(\epsilon_{1}(i,t)+\epsilon_{2}(i,t)+h(i,t))s(t)&\text{if}\
s(t)>0\\\ \epsilon_{3}(i,t)&\text{if}\ s(t)=0.\\\
\end{array}\right.\end{split}$ (1)

$s(t)$ is a latent variable representing the actual number of spots, groups or
composite of the Sun. It cannot be directly observed but its mean may be
estimated based on the observations of the network and later used as a proxy
for $s(t)$. $\epsilon_{1}$ denotes the short-term error representing counting
errors and variable seeing conditions. We assume that
$\mathbb{E}(\epsilon_{1}(i,t))=0$ where $\mathbb{E}$ denotes the expectation
sign. Its variance and tail, that also vary in time, are studied extensively
in Mathieu et al., (2019). $\epsilon_{2}$ represents a long-term error. We are
interested in estimating and monitoring its mean, denoted by $\mu_{2}(i,t)$,
which represents the bias of the stations. At a later stage, we may also study
the variance of $\epsilon_{2}$. We also introduce the random variable
$h(i,t)$, where $h$ is a function of time which varies at a slower pace than
$\epsilon_{2}$. This quantity may be associated to the background level of the
stations (accounting e.g. for differences of instruments or counting
methodologies of the stations) and is not a target of our monitoring
procedures. The errors $\epsilon_{1}$, $\epsilon_{2}$ and $h$ vary on
different time-scales. They affect $s(t)$ in a multiplicative way (Chang and
Oh,, 2012). Finally, $\epsilon_{3}$ captures effects like short-duration
sunspots during solar minima (a detailed study about the distribution of
$\epsilon_{3}$ can be found in Mathieu et al., (2019)).
The random variables $\epsilon_{1}$, $\epsilon_{2}$, $\epsilon_{3}$ and $h$
are assumed to be continuous and $\epsilon_{1}$, $\epsilon_{2}$,
$\epsilon_{3}$, $h$ and $s(t)$ to be jointly independent. Note that
$\epsilon_{1}$, $\epsilon_{2}$ and $\epsilon_{3}$ would be equal to zero and
$h$ be equal to one for a station that would be — in absence of any
measurement errors —perfectly aligned with the solar signal.

###  2.3 Long-term bias

As mentioned earlier, we are not so much interested in the short term error
$\epsilon_{1}$ and our monitoring aims mostly at the component $\epsilon_{2}$
and more specifically its mean $\mu_{2}$. The even longer term levels of the
stations expressed by $h$ are simply estimated and removed from the data. To
this end, we isolate the long-term error from the other components of the
model in the step-wise approach described below. We first divide the
observations by scaling factors to roughly compensate for different observing
conditions: $Z_{i}(t)=\frac{Y_{i}(t)}{\kappa_{i}(t)}$. These piece-wise
constant scaling factors $\kappa_{i}(t)$ are computed as the slope of the
ordinary least-squares regression between the observations of the stations and
the median of the observations ($\underset{1\leq i\leq
N}{\text{med}}Y_{i}(t)$) on periods of 8 months for $N_{s}$, 14 months for
$N_{g}$ and 10 months for $N_{c}$. These values are selected by a statistical-
driven study based on the Kruskal-Wallis test (Kruskal and Wallis,, 1952) that
is completely described in the section 6.3 of Mathieu et al., (2019).
Afterward, we compute $M_{t}$, a robust proxy for $s(t)$ based on the median
of the rescaled observations:

$M_{t}=\underset{1\leq i\leq N}{\text{med}}Z_{i}(t).$ (2)

Motivated from (1), the observations $Y$ are then divided by $M_{t}$ to remove
the main influence of the solar signal. They are also smoothed by a moving-
average (MA) filter, represented by a $\star$ in the following equation. This
smoothing process untangles $eh$, $eh=\epsilon_{2}+h$, from the short-term
error $\epsilon_{1}$:

$\widehat{eh}(i,t)=\left(\frac{Y_{i}(t)}{M_{t}}\right)^{\star}~{}\
\text{when}\ M_{t}>0,$ (3)

where $\widehat{eh}$ denotes the estimator of the mean of $eh$, which is used
as a proxy for $eh$.

Figure 1: Long-term bias, $\hat{\mu}_{2}(i,t)$ for $N_{s}$, in the station
Locarno (Switzerland) over its observing period. The $\hat{\mu}_{2}$s are
smoothed on 27 days (dotted line) to allow the detection of high-frequency
shifts and smoothed on 365 days (plain line) to emphasize the low-frequency
deviations. $M_{t}$ is also represented in the lower plot as an estimation of
the actual number of spots. This figure clearly shows the eleven-year solar
cycle (Hathaway,, 2010) that is intrinsic to the signal.

To analyze the various deviations of the data, different MA-filter window
lengths may be used in (3). The low-frequency shifts such as persisting drifts
are first studied at a yearly scale (i.e. with a window length of 365 days).
Then, a window of length equal to 27 days will also be used to examine the
high-frequency deviations such as sudden jumps. This value of 27 days
corresponds to a physical scale of the data: one solar rotation. It appears to
be sufficiently high to overcome the effects of the short-term regime, as
demonstrated in Mathieu et al., (2019). Finally, the levels of the stations
are separated from the long-term bias by applying once again a MA smoothing
process denoted by $\star\star$:

$\begin{split}\hat{\mu}_{2}(i,t)=\widehat{eh}(i,t)-\widehat{eh}^{\star\star}(i,t),\end{split}$
(4)

where the MA-filter window length should be larger than those of (3). It is
selected here at eleven years or one solar cycle, a physical value that is
larger than the time-scales of the long-term error $\epsilon_{2}$ considered
here. It also seems appropriate since the location of the observatories or
their telescope are unlikely to change much over time.
As an illustration, the $\hat{\mu}_{2}$s smoothed on 27 and 365 days are
represented in Figure 1 for a particular station. Since the $\hat{\mu}_{2}$s
are a (smoothed) ratio, they depend on the level of the solar signal. This is
clearly shown by comparing the lower and upper plots in the figure.

##  3 Ingredients of the method

In this section, the complete monitoring procedure depicted in Figure 2 is
explained. It is intentionally presented in a generic framework to allow the
application of the method on the number of spots $N_{s}$, groups $N_{g}$ as
well as composites $N_{c}$. There are three phases: (I) estimation of the in-
control (IC) parameters of the data, (II) construction and use of the
monitoring procedure and (III) identification of out-of-control patterns.
Phase I contains two steps. At first a subset of stations is selected from the
panel which follows closely the median signal $M_{t}$ mentioned in Section 2.
This pool of stations is then used as a proxy for IC series in the
nomenclature of Qiu and Xiang, (2014). They are used to determine the IC
patterns (mean and variance) of the data and to provide the basis of the
block-bootstrap procedure in Phase II. After standardizing all series by the
IC patterns, the CUSUM control chart is calibrated in phase II by a block
bootstrap procedure from the pool of IC series. The scheme is then applied to
the data for the monitoring. In Phase III, support vector machine (SVM)
procedures predict the shifts size and shape on sub-series detected as out-of-
control by the CUSUM, for easier problem diagnostic.

Selection of IC processes Ia Estimation of IC patterns Ib Design of the chart
and monitoring of standardized data II Estimation of shifts size/shape by SVMs
III Figure 2: Pipeline of the procedure. It is in particular in the dashed
blocks that we contribute with completely new ingredients of the method.

###  3.1 Phase I: Estimation of the IC longitudinal patterns

In this phase I, we automatically construct a subset of IC stations from the
panel and estimate the IC patterns of the data.

####  3.1.1 Phase Ia: Selection of the IC processes

In a first stage, we need a subset (pool) of stations whose observations
follow closely the median signal $M_{t}$. In order to find them, we calculate
a stability criterion on each station. This criterion is based on a robust
version of the mean squared error (MSE) of $\hat{\mu}_{2}$:

$STB(i)=\underset{1\leq t\leq
T}{\text{med}}[\hat{\mu}_{2}(i,t)]^{2}+\underset{1\leq t\leq
T}{\text{iqr}}\hat{\mu}_{2}(i,t),$ (5)

where $\underset{1\leq t\leq T}{\text{iqr}}\hat{\mu}_{2}(i,t)$ and
$\underset{1\leq t\leq T}{\text{med}}\hat{\mu}_{2}(i,t)$ denote respectively
the interquartile range (IQR) and the median of the $\hat{\mu}_{2}(i,t)$ over
the time for a given station $i$. Using these values, we can then cluster the
stations and choose the cluster with the lowest values to form what we call
the pool of IC stations. For this purpose, we use the $k$-means clustering
(Lloyd,, 1957; MacQueen,, 1967) with two clusters.
The pool contains deviations that will be called _disparities_ in the
following, to be distinguished from the deviations that are supposed to be
actually detected by the method. The _disparities_ are expected to be
typically smaller and less frequent than the deviations occurring in the
stations not comprised in the pool. They should be included in the design of
the chart otherwise the scheme would be over-sensitive.

The pool suffers in addition from deviations that are of similar magnitude as
those of the out-of-control (OC) processes. To cope with this and preserve the
detection power of our scheme, we also apply a Shewhart chart (Shewhart,,
1931) with _adaptive_ confidence intervals on the data. We remove the IC
observations that do not fall into one standard deviation around the cross-
sectional mean ($\frac{1}{N}\sum_{i=1}^{N}\hat{\mu}_{2}(i,t)$). Note that we
also test the method with two standard deviations instead of one and the
results were similar. We emphasize that this adaptive Shewhart would not be a
substitute for our control scheme: it only removes the largest deviations at
each time without taking into account the history of the observations.
Therefore, contrarily to our method, it cannot detect the small and persistent
shifts.

####  3.1.2 Phase Ib: Estimation of the mean and the variance of the IC series

We denote by $\mu_{0}(t)$ and $\sigma_{0}^{2}(t)$ respectively the mean and
the variance of the $\hat{\mu}_{2}$ of the pool. Those are estimated by the
empirical mean and variance using nearest neighbours (K-NN) regression method:

$\begin{split}\hat{\mu}_{0}(t)=&\frac{1}{\Delta(t)}\sum_{t^{\prime}=t-\Delta(t)/2}^{t+\Delta(t)/2}\frac{1}{N_{IC}}\sum_{i_{ic}=1}^{N_{IC}}\hat{\mu}_{2}(i_{ic},t^{\prime})\quad
s.t.\quad K=\Delta(t)N_{IC}\\\
\hat{\sigma}_{0}^{2}(t)=&\frac{1}{\Delta(t)}\sum_{t^{\prime}=t-\Delta(t)/2}^{t+\Delta(t)/2}\frac{1}{N_{IC}}\sum_{i_{ic}=1}^{N_{IC}}\left(\hat{\mu}_{2}(i_{ic},t^{\prime})-\hat{\mu}_{0}(t)\right)^{2}\quad
s.t.\quad K=\Delta(t)N_{IC},\end{split}$ (6)

where $i_{ic}$ denotes the index of a station of the pool. With K-NN
regression, the temporal window $\Delta(t)$ can be adjusted to compensate the
missing values of the stations, such that $\hat{\mu}_{0}(t)$ and
$\hat{\sigma}_{0}^{2}(t)$ are always computed on the same number ($K$) of
observations. For appropriate use in the CUSUM chart statistics (to be defined
in (8) below), the data must be standardized by the IC mean and variance (as
in (7) below). Hence the number of nearest neighbors $K$ is selected to obtain
the “best” standardization of the complete panel, in the sense that their
empirical mean becomes close to zero and their empirical variance close to
one. Then, $\Delta(t)$ is chosen in time direction such that
$K=\Delta(t)N_{IC}$.

###  3.2 Phase II: Monitoring

We now turn our attention to monitoring the entire panel. As a reminder, we
are analyzing long-term biases denoted by $\hat{\mu}_{2}$. Using the IC mean
and standard deviation $\hat{\mu}_{0}(t)$ and $\hat{\sigma}_{0}(t)$, we
standardize the (IC and OC) stations to be able to use common monitoring
criteria:

$\hat{\epsilon}_{\hat{\mu}_{2}}(i,t)=\frac{\hat{\mu}_{2}(i,t)-\hat{\mu}_{0}(t)}{\hat{\sigma}_{0}(t)}.$
(7)

Let us now focus on one station (drop the index $i$). We would like to detect
indications of patterns which may relate to problems at the station. This
includes persistent or gradual deviations (shifts or trends) and oscillating
patterns as they may occur when the observatory is used by a rotating pool of
observers each of whom has their own particular way of working. A method for
accumulating small and gradual deviations is to aggregate them over time. A
well known method for doing so in the context of statistical process control
is the cumulative sum (CUSUM) chart (Page,, 1961). The two-sided CUSUM chart
applied on the residuals writes as:

$\begin{split}C_{j}^{+}=max(0,C_{j-1}^{+}+\hat{\epsilon}_{\hat{\mu}_{2}}(t)-k)\\\
C_{j}^{-}=min(0,C_{j-1}^{-}+\hat{\epsilon}_{\hat{\mu}_{2}}(t)+k),\end{split}$
(8)

where $j\geq 1$, $C_{0}^{+}=C_{0}^{-}=0$ and $k>0$ is the allowance parameter
(Qiu,, 2013).
This chart gives an alert if $C_{j}^{+}>h^{+}$ or $C_{j}^{-}<h^{-}$, where
$h^{-}$ and $h^{+}$ are the control limits of the chart. Since the
distribution of the residuals is almost symmetric, we use $h=h^{+}=-h^{-}$.

High deviations may affect the series. Those lead to high values of the CUSUM
statistics which may stay in alert for longer periods than the actual
durations of the shifts. Therefore, in case of too high (resp. too low)
values, we set the chart to a maximal value $2h$ (resp. $-2h$). Hence
$|C^{+}_{j}|,|C^{-}_{j}|\leq 2h$.

####  3.2.1 Design of the chart

As it is clear from the nature of the data, the series to be monitored have a
considerable degree of autocorrelation even when they are in control. We
therefore need a method for determining the control limit of the chart that
takes autocorrelation into account. The block bootstrap (BB) method does this.
It is based on constructing a bootstrap reference distribution by resampling
blocks of data and thereby preserving the autocorrelation of the series.
The control limit ($h$) is adjusted here by a searching algorithm until a pre-
specified rate of false positives (evaluated using the IC average run length,
$ARL_{0}$) is reached with the desired accuracy. The algorithm is explained in
details in Appendix LABEL:app:design (in the supplementary material) and works
as follows. A target shift size, $\delta_{tgt}$, is first selected (it can
also be estimated from the OC series as explained in Appendix
LABEL:app:shift_size). The allowance parameter is specified to
$k=\delta_{tgt}/2$. Then, the actual $ARL_{0}$ is evaluated for an initial
value of the control limit on IC data that are sampled from the pool by the BB
procedure. If the actual $ARL_{0}$ is inferior (resp. superior) to the pre-
specified $ARL_{0}$, the control limit of the chart is then increased (resp.
decreased). This algorithm is iterated until the actual $ARL_{0}$ reaches the
pre-specified $ARL_{0}$ at the desired accuracy.

As theoretically demonstrated in Lahiri, (1999), BB methods using non-
overlapping blocks and random block lengths are more variable than those based
on overlapping blocks and constant lengths. Therefore, we select the popular
moving BB (MBB) (Kunsch,, 1989; Liu and Singh,, 1992) to obtain the best
performances.
Since the BB preserves the serial correlation of the data inside the blocks,
the length of the blocks should be selected appropriately. Large blocks
usually model the autocorrelation of the data properly but at the same time do
not represent well the variance and the mean of the series. And conversely.
Using the method described in Appendix LABEL:app:block_length, the block
length may be selected as the first value such that the MSE of the empirical
autocorrelation of the $\hat{\mu}_{2}$ becomes stable (it corresponds to the
“knee” (or elbow) (Satopaa et al.,, 2011) of the MSE curve). This value
intuitively corresponds to the smallest length which is able to represent the
main part of the autocorrelation of the series.

The data also contain missing values. Among them, there are mainly large gaps
since the smoothing processes in (3) removes the shortest gaps of the series.
As the observing conditions could be different after a large amount of missing
values (different weather conditions or instruments), we restart the scheme
after each gap ($C_{j}^{+}=C_{j}^{-}=0$). Blocks composed only of missing
values are not used to design the chart. This may happen when some stations
contain very few observations on the period studied here, either because they
are ancient and stopped observing at the beginning of the period or because
they have started observing only recently.

###  3.3 Phase III: Estimation of the sizes and shapes of the shifts using
SVMs

The CUSUM gives an alert when a deviation is detected in the data but does not
provide information about the characteristics (shape and size) of the shift.
Such information is however valuable to assign possible causes to the shift or
to adapt the type of alerts that is sent back to the observers. To that end,
Cheng et al., (2011) appended a support vector regression (SVR) to the CUSUM
to predict the magnitude of the shift after each alert. Their method is
however designed for independent and identically normally distributed data
that only experience jumps. Therefore, we extend Cheng et al., (2011) and
design a method that is effective to detect the sizes _and_ the shapes of the
deviations in the sunspot number data. This is achieved by a SVM classifier
(SVC) (Burges,, 1998) in addition to a SVR on top of the chart.

####  3.3.1 Input vector

When an alert is triggered, the $m$ most recent observations of the stations
are fed into the SVR and SVC which then predict the size and the shape of the
deviation at the origin of the alert, as explained in the next subsection. In
particular, the SVR prediction model writes as:

$\hat{\delta}=f(V_{t^{\prime}})=f(\hat{\epsilon}_{\hat{\mu}_{2}}(t^{\prime}-m+1),\hat{\epsilon}_{\hat{\mu}_{2}}(t^{\prime}-m+2),...,\hat{\epsilon}_{\hat{\mu}_{2}}(t^{\prime})),$
(9)

where $t^{\prime}$ denotes the time of the alert and $V_{t^{\prime}}$
represents the input vector, i.e. a sequence containing the last $m$
observations of the series.

The length $m$ of the input vector should thus be sufficiently large to
contain the starting point of most of the deviations while maintaining the
computing efficiency of the method. Large shifts are often quickly detected by
the chart (short OC run length) while the smallest shifts may be identified
only after a certain amount of time (long OC run length). Therefore, the
latter require larger input vectors than the former. Here, $m$ is selected as
an upper quantile of the OC run length distribution for a shift size equal to
$\delta_{tgt}$, as explained in Appendix LABEL:app:input_vector. Hence, $m$
should be sufficiently large to allow the identification of shift sizes that
are superior or equal to $\delta_{tgt}$.
As the SVM procedures do not support missing values, we have to impute them.
Missing observations occurring at the beginning of $V_{t^{\prime}}$ are simply
replaced by the first valid observation encountered, while the “intermediate”
gaps are filled by a linear interpolation. However, when there are too many of
them, the analysis makes no sense. We decide to only analyze input vectors
which have at least 20% of non-missing values.

####  3.3.2 Support vector regression

The support vector machine (SVM) (Vapnik,, 1998) is a supervised machine-
learning procedure, here used as a robust classifier and regressor to predict
the shape and size of the deviations. The method has a strong theoretical
basis that takes root in the optimization theory. It is able to perform
efficiently non-linear classification or regression using a kernel trick that
implicitly maps the data into a high dimension where the non-linear problem
becomes linear. We only introduce the SVR in the following, since the SVC can
be expressed with a similar framework. Smola and Schölkopf, (2004) may also be
consulted for more detailed explanations.
We denote by { $\bm{x}^{j},\delta^{j}|j=1,2,...,M$ } the $M$ training pairs.
$x\in\mathcal{R}^{m}$ represents a training input, i.e. a series of $m$
observations that contains a deviation and $\delta\in\mathcal{R}$ is its
corresponding output, the size of the deviation. The SVR aims at estimating
the continuous regression function relating the deviating observations to the
size of the shift, $f(\bm{x})$, based on the training pairs. This function
writes as:

$f(\bm{x})=\bm{w}^{T}\phi(\bm{x})+b,$ (10)

where $\phi$ is the non-linear function mapping the input data into the high
dimensional feature space, where the regression may be expressed into a
simpler linear problem. The coefficients $\bm{w}$ and $b$ are then estimated
during the training by solving the following optimization problem:

$\underset{\bm{w},b}{\text{min}}\frac{1}{2}||\bm{w}||^{2}+\lambda\frac{1}{M}\sum_{j=1}^{M}L_{\epsilon}(\delta^{j},f(\bm{x}^{j})),$
(11)

where

$L_{\epsilon}(\delta^{j},f(\bm{x}^{j}))=\begin{cases}\begin{tabular}[]{cl}$|\delta^{j}-f(\bm{x}^{j})|-\epsilon$&$\quad|\delta^{j}-f(\bm{x}^{j})|\geq\epsilon$\\\
$0$&otherwise.\end{tabular}\end{cases}$ (12)

In the objective function of (11), the parameter $\lambda$ represents a trade-
off between misclassification and regularization whereas $\epsilon$ in the
loss function of (12) is equivalent to an approximation accuracy, i.e. errors
below $\epsilon$ are neglected. This optimization problem may be rewritten
with Lagrange multipliers as a dual problem and easily solved in the input
space thanks to the introduction of a kernel function
$K(\bm{x^{j}},\bm{x^{k}})=\phi(\bm{x^{j}})\phi(\bm{x^{k}})$. This kernel
function is an important hyper-parameter of the method that should be
carefully selected. After testing different kernels, we choose the radial
basis function, to obtain the best prediction results.

####  3.3.3 Creation of the training and testing sets

The training and testing sets are constructed by simulations since only a
limited amount of observations are available. As the SVM procedures are
supposed to predict the characteristics of the shift after an alert has been
raised by the CUSUM, we generate sets of series that will be first monitored
by the control scheme before reaching the SVMs. When an alert will be
triggered by the CUSUM, the $m$ last values of the series will be assembled
and used as an input vector for the SVMs. Hence, we create series that are
initially longer than $m$. Those are generated from the IC data by BB. To
ensure the efficiency and the generalization of the predictions, we then add
various deviations with different sizes and shapes on top of the series.
_Shift sizes_ The magnitudes of the shifts, $\delta$, are first randomly
sampled from two half-normal distributions (Evans et al.,, 2000) supported by
$[-\infty,...,-\delta_{tgt}]$ and $[\delta_{tgt},...,\infty]$ respectively. We
select the scale parameter of the half-normals equal to 3.5, a value that is
sufficiently high to reproduce the highest values/deviations observed in the
data.
_Shift shapes_ : For each $\delta$, a series $x_{ic}$ of length $T^{\prime}$
is generated from the IC pool by the BB. Here, we choose $T^{\prime}=500$.
Three types of general deviations are then artificially constructed on top of
the series:

  1. 1.

jumps: $x(t)=x_{ic}(t)+\delta$ ;

  2. 2.

drifts with varying power-law functions:
$x(t)=x_{ic}(t)+\frac{\delta}{T^{\prime}}(t)^{a}$, where $a$ is randomly
selected in the range $[1.5,2]$ ;

  3. 3.

oscillating shifts with different frequencies: $x(t)=x_{ic}(t)\sin{(\eta\pi
t)}\delta$, where $\eta$ is randomly selected in the range
$[\frac{\pi}{m},\frac{3\pi}{m}]$.

These deviations are not data-specific to guarantee the generalization of the
predictions of the characteristics of the shifts.
_Time of the shift_ : In the data, the shifts may happen not immediately but
after an initial IC period. Therefore, we also start the monitoring after a
random delay in the range $[m,3m/2]$, to train the methods at identifying
shifts appearing anywhere within the input vector. Note that the SVMs as well
as the control chart should be started after $m$ observations are gathered.

The SVR is trained on these constructed sets to predict the size of the
deviations in the continuous range
$[-\infty,...,-\delta_{tgt},\delta_{tgt},...,\infty]$. In practice, we observe
that the SVR generalizes well and can make predictions on $\mathbb{R}$ even if
it was only trained in a smaller range of interest. The SVC also learns on the
same sets to identify three different shapes: jumps, drifts and oscillating
shifts. If a wide range of deviations are simulated, only three classes are
therefore involved in the classification problem.

##  4 Monitoring the composite sunspot index Nc

In this section, we use our methodology to solve the monitoring problem of the
sunspot numbers. We do so for the composite $N_{c}=N_{s}+10N_{g}$ (the same
approach also works for the two components, $N_{s}$ and $N_{g}$ and is
presented in the supplemental material). We first study the low-frequency
deviations on data that have been smoothed on a year to extract and analyze
low-frequency patterns such as trends and persistent shifts. Then, we examine
on data that have been smoothed on 27 days (one solar rotation) the higher
frequency patterns such as sudden jumps. The section ends with an example of a
multi-scale monitoring applied to a particular observatory.

###  4.1 Lower frequency monitoring

In the first step of low-frequency monitoring of $N_{c}$, we smooth the long-
term bias ($\hat{\mu}_{2}$) with a window length of one year as described in
Section 2.3. As explained in Section 3.1.1, the network of stations is first
reduced to a pool of 119 in-control (IC) stations. In the next step, we
extract the IC mean and standard deviation using the K-NN regression described
in Section 3.1.2. The selection mechanism finds $K=4600$ for this step. The
resulting mean and standard deviation are then used to standardize all series.

In the second stage, we use the block bootstrap method described in Section
3.2.1 to calibrate the CUSUM chart at an average run length of 200. This
requires choosing the block length first. In our situation, a choice of two
solar rotations (54) appears appropriate. It is longer than the lifetime of
most sunspots but not too long for practical use. The calibration then leads
to a control limit of $h=19$ and a target shift size of $\delta_{tgt}=1.5$.

(a)

(b)

Figure 3:  (a) Upper panel: the residuals
$\hat{\epsilon}_{\hat{\mu}_{2}(i,t)}$ for $N_{c}$ smoothed on one year from
the KZ over the period studied (1981-2019). In addition to their disparities,
the residuals also contain the actual deviations of the station, which have
been removed for the design of the chart as explained in Section 3.1.1. Middle
panel: the (two-sided) CUSUM chart statistics applied on the residuals in
_square-root scale_. The control limits of the chart are represented by the
two horizontal thick lines. Lower panel: the characteristics of the deviations
predicted by the SVR and SVC after each alert. (b) Similar figure for MT over
the same period.

Finally, the support vector method for extracting and classifying out-of-
control patterns is deployed. It is composed of a SVR to predict the size of
the shifts and a SVC to classify the shape of the encountered deviations. We
obtain them by creating a set of artificial series of 500 values generated
from the IC pool by the BB. These give us series as we would observe in
reality including correlations. We then artificially add jumps, trends, and
oscillating shifts to them as described in Section 3.3.3. These series are
then fed to the CUSUM chart which identifies the out-of-control observations.
When an alert is triggered by the chart, an input vector containing the $m$
last observations of the series is assembled. This input vector is then
analyzed by the SVR and SVC for predicting the characteristics of the shift.
In our case, we harvested 63000 such series from the IC pool and we enriched
them with artificial patterns. We then calibrated SVR and SVC models by
splitting this set into a training set (80%) and a testing set (20%).

The length of the input vector is specified here at $m=80$, as explained in
Section 3.3.1. This value of 80 corresponds to the 90-th quantile of the OC
run length distribution, for a shift of size $\delta_{tgt}$. The other
parameters of the support vector machines are automatically selected from a
searching interval to obtain the best prediction results. Those are evaluated
using the mean absolute percentage error (MAPE) for the regression and the
accuracy for the classification problem, see Appendix LABEL:app:performances.
With this method, the regularization parameter $\lambda$ of the classifier and
regressor is set to 13 and the accuracy error $\epsilon$ of the SVR is fixed
at 0.001. The performances of the SVMS are presented in Appendix
LABEL:app:perfo_drifts. Overall, they are sufficient to achieve our goals:
identify the origins of the deviations.

Figure 3 shows results for two stations labeled KZ333The Kanzelhöhe
Observatory in Austria. and MT444The National Observatory of Japan, in Mitaka
(Tokyo). The observatory of KZ is rather stable and belongs to the IC pool. It
relies on a stable team of well trained observers. The observatory MT is
generally less stable and shows a severe downward trend around 1998. The cause
of this downward trend could be traced to the replacement of visual counts
from the direct optical solar image by automatic computerized counts based on
digital images from a CCD camera. Given the image sensor technology then
available, the spatial resolution of the images was limited, and many small
spots that were fully detected in earlier visual observations were not
detected anymore by the new equipment.

###  4.2 Higher frequency monitoring

The method described above can also be applied to the biases ($\hat{\mu}_{2}$)
smoothed on a shorter time window such as the duration of a solar cycle (27
days). Here the selection of the IC pool yields 100 stations and the number of
nearest neighbors comes out to $K=2400$. This number is smaller than before
since we are working at a higher frequency. The block bootstrap and SVMs with
the same settings as above can be used to calibrate the CUSUM chart. For
$\delta_{tgt}=1.4$, the control limit of the chart is selected at $h=13$ to
obtain an average run length of 200. The length of the input vector is fixed
here at $m=70$, which corresponds to the 90-th quantile of the OC run length
distribution.

(a)

(b)

Figure 4:  (a) Upper panel: the residuals
$\hat{\epsilon}_{\hat{\mu}_{2}(i,t)}$ for $N_{c}$ smoothed on 27 days for KO
over the period studied (1981-2019). In addition to their disparities, the
residuals also contain the actual deviations of the station, which have been
removed for the design of the chart as explained in Section 3.1.1. Middle
panel: the (two-sided) CUSUM chart statistics applied on the residuals in
_square-root scale_. The control limits of the chart are represented by the
two horizontal thick lines. Lower panel: the characteristics of the deviations
predicted by the SVR and SVC after each alert. (b) Similar figure for SM over
the same period.

Figure 4 shows the methodology applied to KO555Name of the observer known to
the authors, kept for privacy. and SM666The observatory of San Miguel in
Argentina.. KO was a Japanese observatory run by a single dedicated observer
whose records (which stopped during 1996) were very stable. On the contrary,
SM is a severely OC station that experiences large known deviations (Mathieu
et al.,, 2019, Figure 12). The large variations observed in SM are likely
caused by the rotation of several observers involved in the counting process.
In some countries, the public observatories have also an educational function.
Their team of regular observers are usually small and are often completed by
student or amateur astronomers that are frequently replaced, which causes
large variations. Unfortunately, we could not identify more precisely the
origin of the deviations since their observations stopped several years ago
and we have not succeeded in contacting them yet. The lack of information is a
common problem we face when investigating past deviations in stations that are
now inactive and is therefore worth mentioning.
As we see in the figure, the biases vary a lot at 27 days, a scale which is
close to the short-term regime. The actual monitoring should be based on a
larger scale otherwise some stations such as SM would receive almost constant
alerts. If a particular deviation is detected at higher scale (such as one
year), it might be interesting however to analyze it at 27 days, to better
identify its origin.

###  4.3 Multi-scale monitoring

(a)

(b)

Figure 5:  (a) Upper panel: the residuals
$\hat{\epsilon}_{\hat{\mu}_{2}(i,t)}$ for $N_{c}$ smoothed on 365 days from FU
over the period studied (1981-2019). In addition to their disparities, the
residuals also contain the actual deviations of the station, which have been
removed for the design of the chart as explained in Section 3.1.1. Middle
panel: the (two-sided) CUSUM chart statistics applied on the residuals in
_square-root scale_. The control limits of the chart are represented by the
two horizontal thick lines. Lower panel: the characteristics of the deviations
predicted by the SVR and SVC after each alert. (b) Similar figure for the
values of $\hat{\epsilon}_{\hat{\mu}_{2}(i,t)}$ smoothed on 27 days in FU.

Figures 3 and 4 display instances of a stable IC station included in the pool
and a typical out-of-control observatory for the high- and low- frequency
monitoring respectively. To better grasp the motivations of a multi-scale
monitoring, the method is applied to the data smoothed on 27 days and one year
of FU777Name of the observer known to the authors, kept for privacy. in Figure
5. The FU station is composed of a single dedicated observer in Japan, who has
observed without interruption since 1968 until today, producing one of the
longest individual series. His observations are included in the IC pool but
yet suffer from recent deviations. In particular, the upward deviation (which
looks like a spike) reported in 2007 in FU (Clette,, 2013) as well as the
downward drift occurring after 2014 are well identified in Figure 5(a). Figure
6 shows a zoom of Figure 5(a) on the time period from 2007 to 2008. After a
progressive upward shift, the station experiences a rapid downward trend over
five days. This trend, which looks like a jump in the whole period view, it
thus correctly classified by the SVC. Even by taking a closer look on the
figure however, it remains difficult to precisely identify the origin of the
shifts on data that are smoothed on a year. By looking at a smaller scale of
27 days in Figure 5(b), we can better characterize the shift in 2007 as a
short event and pinpoint its location. After investigations, this deviation
appears to be related to a small over-count that appeared in early 2007 (three
groups were reported in FU while most of the network only observed two groups)
while the drift might be associated to the health condition of the observer.
Note that the long-term biases are not defined (i.e. set to missing values)
when the median of the network is equal to zero, see (3). This regime
corresponds to those of the variability at minima, represented by
$\epsilon_{3}$. Due to the smoothing procedure of (3), the deviations that
appeared close to solar minima, such as the jump in FU, are thus particularly
visible.
As shown in the figures, the monitoring and the SVM procedures can cope with a
large variety of shifts ranging from small and persistent deviations to large
oscillating shifts. The procedures automatically detect major deviations
recently discovered by hand as mentioned above. More identified prominent
deviations as well as results for other stations are shown in Appendix
LABEL:app:figures. In addition, the chart also unravels many other shifts,
typically smaller, that are otherwise difficult to identify.

Figure 6: Sizes and shapes of the deviations (taken from Figure 5a) predicted
by SVMs in FU over 2007-2008.

Note that the figures represent past observations, which have not been
monitored by any control scheme. Consequently, the stations may stay in alert
for long consecutive periods. If this method is applied on future
observations, any major deviation will be promptly corrected. Therefore we
expect better results for data that have already been monitored.

##  5 Conclusion and perspectives

We presented a nonparametric control scheme to monitor challenging and
important datasets related to the observations of sunspots across a wide
network of stations. The approach allows us to deal with the missing values,
autocorrelations, and non-normality of the data, to detect and classify
station-related anomalies on different time scales. The procedure is based on
a particular choice of methods for smoothing, robust anomaly detection, and
anomaly classification. Other methods exist for these steps; yet we believe
that our choices are particularly suitable for the problem at hand.

The features of our approach are multi-scale smoothing, CUSUM charting, SVM
classification and detailed graphical displays. They also include an automatic
pre-selection of an in-control pool and the powerful calibration of the chart
using block bootstrap procedures. The associated advantages are robustness,
flexibility, automation, and guided interpretation of results. The method
allows us to detect and identify the causes of major deviations that occurred
in the series. We have seen that monitoring on at least two-time scales is
essential to capture these anomalies. Some patterns first attract interest on
a long-term scale but it is at the short-term scale that their potential root-
causes can be suggested. The method also identifies a wide range of deviations
unseen in previous analyses. Most of them have not been related to specific
causes yet but will soon be investigated.

This automated method allows us and the researchers at the Royal Observatory
who are in charge of producing the International Sunspot Number to have a
harmonized view across the network of stations. It provides a way to give
specific and targeted advice to the observers. As demonstrated in this paper,
the method also delivers easy to interpret graphical displays which facilitate
root cause analysis of deviations. The complete re-examination of past data of
the whole panel has just started. When they will be finished, these analyses
will allow us to arrive at a cleaner data stream and to release an improved
version of the International Sunspot Number. Additionally, the implementation
of the method in the continuous surveillance of future observations will lead
to a faster detection and identification of inconsistencies, their elimination
by better observer training or equipment maintenance, and finally to a more
precise determination of the sunspot numbers in the future.

Our methodology has the potential to adapt to other type of panel data such as
those observed in the manufacturing or financial industry, which remains to be
investigated.

SUPPLEMENTARY MATERIAL

Python package (codes)

The subset of data and the codes that we used in this paper are available at
https://github.com/sophiano/SunSpot.

Algorithms

The pseudo-algorithms to design the CUSUM chart, to select the target shift
size, to choose the block length and the length of the input vector are
explained in the supplementary material in Appendix A.

Performance criteria

The performances of the support vector procedures for the high and low
frequency monitoring are displayed in the supplementary material in Appendix
B.

Additional figures

Additional analyzes and figures are also provided in the supplementary
material in Appendix C.

##  6 Acknowledgments

The first author gratefully acknowledges funding from the Belgian Federal
Science Policy Office (BELSPO) through the BRAIN VAL-U-SUN project
(BR/165/A3/VAL-U-SUN). S. Mathieu, L.Lefèvre and F.Clette also acknowledge
financial support from the International Space Science Institute (ISSI, Bern,
Switzerland) via the International Team 417 “Recalibration of the Sunspot
Number Series”, chaired by M. Owens and F.
Clette888https://www.issibern.ch/teams/sunspotnoser/. Computational resources
have been provided by the Consortium des Équipements de Calcul Intensif
(CECI), funded by the Fonds de la Recherche Scientifique de Belgique
(F.R.S.-FNRS) under Grant No. 2.5020.11 and by the Walloon Region. The authors
also thank Hisashi Hayakawa and all observers that provided their help in the
search for the origins of the deviations.

## References

  * Burges, (1998) Burges, C. (1998).  A tutorial on support vector machines for pattern recognition.  Data Mining and Knowledge Discovery, 2:121–267.
  * Chang and Oh, (2012) Chang, H.-Y. and Oh, S.-J. (2012).  Does correction factor vary with solar cycle?  Journal of Astronomy and Space Sciences, 29(2):97–101.
  * Cheng et al., (2011) Cheng, C.-S., Chen, P.-W., and Huang, K.-K. (2011).  Estimating the shift size in the process mean with support vector regression and neural network.  Expert Systems with Applications, 38(8):10624–10630.
  * Clette, (2013) Clette, F. (2013).  Private Communication: Talk presented at 3rd Sunspot Number Workshop (Tucson, USA).
  * Dudok de Wit et al., (2016) Dudok de Wit, T., Lefèvre, L., and Clette, F. (2016).  Uncertainties in the sunspot numbers: estimation and implications.  Solar Physics, 291(9-10):2709–2731.
  * Ermolli et al., (2013) Ermolli, K., Matthes, K., Dudok de Wit, T., Krivova, N., Tourpali, K., Weber, M., Unruh, Y., Gray, L., Langematz, U., Pilewskie, P., Rozanov, E., Schmutz, W., Shapiro, A., Solanki, S., and Woods, T. (2013).  Recent variability of the solar spectral irrandiance and its impact on climate modelling.  EGU publication : Atmospheric Chemistry and Physics, 13:3945–3977.
  * Evans et al., (2000) Evans, M., Hastings, N., and Peacock, B. (2000).  Statistical Distributions.  Wiley, New-York, 3rd edition.
  * Haigh, (2002) Haigh, J. (2002).  The effects of solar variability on the Earth’s climate.  Philosophical Transactions of the Royal Society: Mathematical, Physical and Engineering Sciences, 361(1802):95–111.
  * Hathaway, (2010) Hathaway, D. H. (2010).  The solar cycle.  Living Reviews in Solar Physics, 7:1.
  * Kruskal and Wallis, (1952) Kruskal, W. and Wallis, W. (1952).  Use of ranks in one-criterion variance analysis.  Journal of the American Statistical Association, 47(260):583–621.
  * Kunsch, (1989) Kunsch, H. R. (1989).  The jackknife and the bootstrap for general stationary observations.  The Annals of Statistics, 17(3):1217–1241.
  * Lahiri, (1999) Lahiri, S. N. (1999).  Theoretical comparisons of block bootstrap methods.  The Annals of Statistics, 27(1):386–404.
  * Liu and Singh, (1992) Liu, R. Y. and Singh, K. (1992).  Moving blocks jackknife and bootstrap capture weak dependence.  In Lepage, R. and Billard, L., editors, Exploring the Limits of Bootstrap, pages 225–248. Wiley, New-York.
  * Lloyd, (1957) Lloyd, S. (1957).  Least squares quantization in PCM.  Technical Report RR-5497 5497, Bell Lab.
  * MacQueen, (1967) MacQueen, J. B. (1967).  Some methods for classification and analysis of multivariate observations.  In Le Cam, L. and Neyman, J., editors, Proceedings of the fifth Berkeley symposium on mathematical statistics and probability, pages 281–297. University of California Press, California, United States.
  * Mathieu et al., (2019) Mathieu, S., von Sachs, R., Delouille, V., Lefevre, L., and Ritter, C. (2019).  Uncertainty quantification in sunspot counts.  The Astrophysical Journal, 886(1):7.
  * Morfill et al., (1991) Morfill, G., Scheingraber, H., Voges, W., and Sonett, C. (1991).  Sunspot number variations -stochastic or chaotic.  In Sonett, C., Giampapa, M., and Matthews, editors, The Sun in Time, pages 30–58. University of Arizona Press, Tucson, United States.
  * Page, (1961) Page, E. S. (1961).  Cumulative sum charts.  Technometrics, 3(1):1–9.
  * Qiu, (2013) Qiu, P. (2013).  Introduction to Statistical Process Control.  CRC Press, Taylor and Francis Inc, 1st edition.
  * Qiu and Xiang, (2014) Qiu, P. and Xiang, D. (2014).  Univariate dynamic screening system: an approach for identifying individuals with irregular longitudinal behaviour.  Technometrics, 56(2):248–260.
  * Satopaa et al., (2011) Satopaa, V., Albrecht, J., Irwin, D., and Raghavan, B. (2011).  Finding a ”kneedle” in a haystack: Detecting knee points in system behavior.  In 2011 31st International Conference on Distributed Computing Systems Workshops, pages 166–171. IEEE.
  * Shewhart, (1931) Shewhart, W. (1931).  Economic Control of Quality of Manufactured Product.  Van Nostrand, 1st edition.
  * Smola and Schölkopf, (2004) Smola, A. and Schölkopf, B. (2004).  A tutorial on support vector regression.  Statistics and Computing, 14(3):199–222.
  * Temmer et al., (2001) Temmer, M., Veronig, A., Hanslmeier, A., Otruba, W., and Messerotti, M. (2001).  Statistical analysis of solar H$\alpha$ flares.  Astronomy and Astrophysics, 375(3):1049–1061.
  * Vapnik, (1998) Vapnik, V. (1998).  Statistical Learning Theory.  Wiley, New-York, 1st edition.
  * Vigouroux and Delache, (1994) Vigouroux, A. and Delache, P. (1994).  Sunspot numbers uncertainties and parametric representations of solar activity variations.  Solar Physics, 152(1):267–274.
  * Wang and Colaninno, (2014) Wang, Y.-M. and Colaninno, R. (2014).  Is Solar Cycle 24 producing more coronal mass ejections than Cycle 23?  The Astrophysical Journal Letters, 784(L27):1–7.


# Shared Loss between Generators of GANs

Xin Wang
Independent Researcher
<EMAIL_ADDRESS>
###### Abstract

Generative adversarial networks are generative models that are capable of
replicating the implicit probability distribution of the input data with high
accuracy. Traditionally, GANs consist of a Generator and a Discriminator which
interact with each other to produce highly realistic artificial data.
Traditional GANs fall prey to the mode collapse problem, which means that they
are unable to generate the different variations of data present in the input
dataset. Recently, multiple generators have been used to produce more
realistic output by mitigating the mode collapse problem. We use this multiple
generator framework. The novelty in this paper lies in making the generators
compete against each other while interacting with the discriminator
simultaneously. We show that this causes a dramatic reduction in the training
time for GANs without affecting its performance.

##  1 Introduction

Generative adversarial networks (GANs) [1] are highly successful in
replicating the implicit probability distribution of any dataset. It has been
widely used in a variety of applications involving generation of realistic
images [2], network analysis [3, 4, 5] and generation of privacy preserving
synthetic datasets [6, 7, 8]. The primary reason for its widespread adoption
is that it produces much more realistic data than its predecessors like
variational autoencoders [9].

Unlike the traditional maximum likelihood estimation based generative methods,
GANs are generative functions which rely on game theory. When provided with a
sample, the discriminator network D(x) tries to determine whether it is drawn
from the generator distribution $p_{model}(x)$ or the real data distribution
$p_{data}(x)$. The generator network G(z;$\theta^{(G)}$) produces artificial
data and tries to deceive D(x) into accepting it as a real data sample. The
confidence with which D(x) classifies it as a real data sample is fed back to
the generator as its training signal. This process continues until the game
reaches a Nash equilibrium, that is $D(x)$ cannot correctly distinguish
between $p_{model}(x)$ and $p_{data}(x)$.

The GANs are trained using gradient descent methods which do not provide any
theoretical guarantee about its convergence [10]. Methods have been proposed
to train GANs more efficiently [11] but there is still no guarantee that the
training will converge $p_{model}(x)$ to $p_{data}(x)$ eventually. One of the
major drawbacks of GANs is the problem of ‘mode collapse’ [12] [13] [14] [15]
[11]. Theoretically, the generator is expected to learn the true data
distribution at convergence. However, in practice, the generator cannot learn
distributions with multiple modes correctly. For data distributions having
multiple modes at $x=M_{1},M_{2},...,M_{k}$, it is observed that
$p_{model}(x)$ has a single mode at $x=M_{i_{1}}$, such that
$i_{1}\epsilon\\{1,2,...,k\\}$.

This problem can be addressed by improving the training methodologies of GANs
[12] [15] [11]. Another method of addressing this problem is to enforce the
GAN to capture the diverse modes through various architectural modifications
[16] [14] [13].

In this work, we propose a new architecture which trains much faster and
converges much better than the existing methods. This architecture is
motivated by the target to achieve better training performance and convergence
in the new architectures for tackling mode-collapse.

This paper is further organized as follows: Section 2 provides a background
about GANs and their working dynamics. In Section 3, we describe our
motivations for pursuing this idea and its merits. Section 4 presents the
methodology and detailed implementation of our idea. Section 5 includes the
experimental results, followed by discussions on the merits and possible
drawbacks of this scheme in Section 6. We conclude this paper in Section 7 and
also suggest some directions of future work.

##  2 Background

Given a generator $G$ and a discriminator $D$, the training process of a GAN
can be formulated as the following min-max game:

$\min_{G}\max_{D}\operatorname{\mathbb{E}}_{x\leftarrow{P_{data}}}[logD(x)]+\operatorname{\mathbb{E}}_{z\leftarrow{P_{z}}}[1-logD(G(z)]$

where $x$ is drawn from data distribution $P_{data}$ and $z$ is drawn from a
prior distribution $P_{z}$. The mapping $G(z)$ induces a generator
distribution $P_{model}$ in data space. GAN alternatively optimizes $D$ and
$G$ using stochastic gradient-based learning. $G$ is therefore incentivized to
map every $z$ to a single $x$ that is most likely to be classified as true
data, leading to the mode collapsing problem. Another commonly asserted cause
of generating less diverse samples in GAN is that, at the optimal point of
$D$, minimizing $G$ is equivalent to minimizing the Jensen-Shannon divergence
(JSD) between the data and model distributions, which has been empirically
proven to prefer to generate samples around only a few modes whilst ignoring
other modes [17] [18]. The working methodology of a GAN is depicted in Fig.1.

Figure 1: Working of a GAN [19]

##  3 Motivation

The challenge of mode collapse faced while training GANs has been tackled by
using multiple generators. Various multi-generator GAN architectures [20] [21]
[22] demonstrate much better coverage of modes of the data distribution using
multiple generators. It has also been theoretically proven that multiple
generators are indeed successful in tackling the mode-collapse problem [23].
When there are multiple generators, the discriminator not only has to guess if
the sample originates from the real data distribution or the generator, but
also has to determine from which generator the sample has originated. So, the
discriminator forces the generators to produce different outputs from the real
data distribution, thus ensuring diversity and evading mode collapse. Multiple
generators are also better at producing better artificial data than GANs with
a single generator [24].

Other multi-agent GAN architectures also include the GMAN [16] architecture
which involves multiple discriminators. It has been shown that using multiple
discriminators can effectively train GANs in a fraction of the time required
originally. However, GANs with multiple discriminators cannot address the
mode-collapse problem.

In this work, we propose a new class of GANs which consists of multiple
generators. These generators not only learn from the training signal provided
by the discriminator but also by competing against each other. We call these
class of GANs Racing-GAN and show that Racing-GANs train much faster than
existing GAN frameworks while tackling mode collapse at the same time.
Previously, the multiple generators collaborated with each other by the means
of shared weights [16] or by having some auxiliary information.
Parallelization of GANs has also been proposed in which multiple generators
are connected to multiple discriminators without any interconnection between
the generators [25]. The novelty of our work lies in initiating a competition
of the generators among themselves to train at a much faster rate.

##  4 Methodology

Racing-GAN involves a competition of multiple generators with each other. This
involves an intricate interconnection of the generators among themselves. In a
model with $k$ generators, $2^{k}$ such interconnections are possible.
However, in this introductory paper we explore the simplest architecture
possible. We connect every generator to a single discriminator. Every
generator generates its own output as shown in Fig. 2. In our experiments, we
consider the output of any one of the generators as the final output. The
choice of this generator can be random, as we show that all the generators
perform equally well at the end of training. All these outputs can also be
combined to produce a single output as in [23]. This may include a voter which
receives the outputs of all the generators and outputs the one it considers to
be the best.

Figure 2: Racing-GAN model with $k$ generators and a discriminator

###  4.1 Detailed Implementation

Previously, multiple methods have been employed to make the outputs dependent
on multiple generators. Some schemes utilize weight sharing [22] while others
use auxiliary information [24]. In Racing-GAN, we modify the loss functions of
the generators to make them compete with each other.

In a model with one Discriminator and $k$ generators, the objective function
of the Discriminator is:

$\max_{D}\operatorname{\mathbb{E}}_{x\leftarrow{P_{data}}}[logD(x)]+\operatorname{\mathbb{E}}_{z\leftarrow{P_{z}}}[\sum_{n=1}^{k}\\{1-logD(G_{n}(z)\\}]$

This loss function forces the Discriminator to not only guess if the input
sample belongs to the real data or not, but also to determine from which
generator an artificial data sample has originated. This causes the
discriminator to influence the generators to produce non-similar images, hence
tackling the mode-collapse problem. The analysis of this loss function is
explained in greater depth in Section 4.2.

The objective function of a generator $G_{i}$ in a model in which all the
generators compete with each other, i.e., a fully connected model is

$\min_{G_{i}}\operatorname{\mathbb{E}}_{z\leftarrow{P_{z}}}[1-log(D(G_{i}(z)+\sum_{j\epsilon\\{1,..,k\\},j\neq
i}\\{max(0,D(G_{i}(z))-D(G_{j}(z)))\\}]$

In this paper, we consider a model in which a Generator $G_{i}$ is connected
to only one Generator $G_{j}$. The objective function of Generator $G_{i}$
coupled with a single Generator $G_{j}\forall i,j\epsilon\\{1,2,...,k\\},i\neq
j$:

$\min_{G_{i}}\operatorname{\mathbb{E}}_{z\leftarrow{P_{z}}}[1-log(D(G_{i}(z)+max(0,D(G_{i}(z))-D(G_{j}(z)))]$

###  4.2 Analysis

The only condition for a GAN to be functional is that the loss functions
should be differentiable. However, we introduce the maximum function in the
loss function of the Generator but it is observed that the model works better
than the original. This phenomenon can be viewed as the introduction of noise
in the loss function of the GAN. Introduction of noise is a common practice in
deep learning [26] to improve learning and generalization. The objective
function of the Racing-GAN can be viewed as variant of addition of noise in
certain iterations. Addition of noise improves learning with gradient descent
and GANs are trained using gradient descent. This intuitively explains the
reason why Racing-GANs have a better training performance than the other
variants of GANs. Since addition of noise during gradient descent also
improves generalization, the loss functions of Racing-GANs face much lower
oscillations after convergence than their counterparts.

All the generators start with random samples from the same probability
distribution. It is expected that their performances will be nearly equal. But
whenever generator $G_{i}$ performs better than generator $G_{j}$, $G_{j}$ is
penalized by the addition of a positive term in its loss function. This causes
$G_{j}$ to update its weights by a greater value, hence causing it to learn
faster.

An interesting observation suggests that, due to the competitive nature of the
generators in Racing-GAN, the loss functions of all the generators are closely
tied to each other and tend to follow the same curve after the first few
iterations. This is beneficial since it causes minimal oscillations in the
loss functions after convergence.

##  5 Experiments

We present the experimental results in two parts. Firstly, we demonstrate that
the Racing-GAN indeed produces good artificial data which are comparable to
other GANs. Then we show that these data are produced after much less training
as compared to other variants of GANs.

###  5.1 Experimental Setup

We have used a simple GAN model for our initial experiments. The generator
objective is to generate a two-dimensional quadratic curve which lies between
two predetermined quadratic curves. For simplicity, we use Racing-GANs with
only two generators, namely $G_{1}$ and $G_{2}$.

For our experiments, we use four frameworks for comparing and highlighting the
efficiency of the induced competition between the generators. These frameworks
are:

  * •

GAN 1: An ordinary GAN with only one Generator $G$

  * •

GAN 2: A GAN with two uncoupled independent Generators $G_{1}$ and $G_{2}$

  * •

GAN 3: A GAN with one Generator $G_{2}$ competing with $G_{1}$ whereas $G_{1}$
is independent of $G_{2}$

  * •

GAN 4 (Racing GAN): Both the Generators $G_{1}$ and $G_{2}$ compete with each
other.

###  5.2 Results

In this section, we present the outputs of the four variants of the GANs
mentioned above after the GANs have reached their respective equilibria. They
are successful in generating quadratic curves between the specified
boundaries. This can be observed in the output of the iteration 8000 of the
GANs in Fig.3 and Fig.4. From these figures, we get a rough idea that multi-
generator frameworks are better than the ordinary GAN framework. This aspect
is demonstrated by analyzing the convergence rate of the loss functions of
these different frameworks in a later section. There we show that GAN 4 is the
best among the variants considered over here.

Figure 3: Generator outputs at iterations 1, 2500, 8000 for an ordinary GAN with single generator  Figure 4: Multiple generator outputs st iterations 1, 2500, 8000 for (from top) GAN 2, GAN 3, GAN 4 Table 1: Iterations required by Discriminator $D$ and Generators $G1$ and $G2$ to converge GAN type | D | G | G1
---|---|---|---
GAN 1 | 9348 | 8107 | –
GAN 2 | 8547 | 4484 | 5872
GAN 3 | 8969 | 5364 | 6910
GAN 4 | 7728 | 6049 | 4944
Table 2: Improvement in training efficiency over the ordinary GAN (GAN 1)  GAN type | D | G
---|---|---
GAN 2 | 8.57% | 44.69%
GAN 3 | 4.05% | 33.83%
GAN 4 | 17.33% | 39.02%

###  5.3 Loss functions

In this section, we plot the loss functions of the Discriminator and the
Generators for each of the four variants of the GANs mentioned above. From
Tables 1 and 2, we observe that the loss functions for the Racing-GAN converge
much faster than the other variants. This shows that training of a Racing-GAN
occurs faster than the other variants. (We consider that equilibrium is
attained when the value of the loss function oscillates within $\pm 1\%$ of
its correct value.)

(a) GAN 1

(b) GAN 2

(c) GAN 3

(d) GAN 4

Figure 5: Loss functions of Discriminator (red), Generator $G_{1}$ (green) and
Generator $G_{2}$ (blue)

While using the GAN to generate synthetic data, the output from any one of the
generators is taken. If a multi-generator GAN is used, the generator producing
the best output will be utilized. Hence, to compare fairly with other GAN
variants we have used the generator which performs the best in case of
multiple generators.

###  5.4 Observations

The observations stated here have been observed over several experimental
iterations over the aforementioned GAN variants.

  * •

From Table 2, it can be observed that a GAN with one competing generator and
another independent generator (GAN 3) performs worse than a GAN with two
independent generators (GAN 2).

  * •

In Figures 4 (a),(b),(c), it is observed that the loss functions of $G_{1}$
and $G_{2}$ closely imitate each other. This is due to the competing nature of
the generators. Both of them want to be as good as the other, leading to a
close overlap of their loss functions.

  * •

The generator of GAN 2 performs slightly better than that of GAN 4 but the
discriminator of GAN 4 converges at least twice as fast as that of GAN 2. This
has been observed in every experimental trial conducted. Hence, we can
conclude that GAN 4 which is a Racing-GAN with two competing generators
perform significantly better than the other variants.

##  6 Discussions

The Racing-GAN is an interesting architecture of generative adversarial
networks which train much faster than its counterparts. It is expected to
solve the mode collapse problem. We expect to see an increased interest in the
research community where various generators are combined in different fashions
to get better results. The Racing-GAN is an approximately differentiable
model, which is inspired from nature, seems to work better than the other
variants. However, one of the basic requirements for GANs to operate is that
the loss function has to be differentiable. This is a contradiction which is
unsolved as of now. An interpretation of the success of Racing-Gan in spite of
being non-differentiable is that the extra term being added to the generator
loss function is being treated as noise. It has been shown that addition of
noise improves generalization and avoids overfitting in deep neural
networks[26].

##  7 Conclusion and Future Work

We present a new model of GAN called a Racing-GAN which trains much faster and
effectively than the existing variants of GANs. This model has multiple
generators which compete with each other besides learning from the training
signal provided by the discriminator. The work presented in this paper is of
practical and experimental nature. We hope to develop a deeper theoretical
understanding of its dynamics in the near future. Other future work might also
include investigating various couplings between the generators and modifying
the loss functions for better results.

## References

  * [1] I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. C. Courville, and Y. Bengio, “Generative adversarial nets,” in _NIPS 2014_.
  * [2] H. Zhang, T. Xu, and H. Li, “Stackgan: Text to photo-realistic image synthesis with stacked generative adversarial networks,” in _Proceedings of the IEEE, Oct 2017_.
  * [3] T. Saha, N. Aaraj, and N. K. Jha, “Machine learning assisted security analysis of 5g-network-connected systems,” _IEEE Transactions on Emerging Topics in Computing_ , 2022.
  * [4] ——, “System and method for security in internet-of-things and cyber-physical systems based on machine learning,” Jun. 23 2022, uS Patent App. 17/603,453.
  * [5] T. Saha and V. Sehwag, “Tv-puf: A fast lightweight aging-resistant threshold voltage puf,” _Cryptology ePrint Archive_ , 2016.
  * [6] I. J. Goodfellow, “Differentially private releasing via deep generative model (technical report),” 2018.
  * [7] T. Saha, N. Aaraj, N. Ajjarapu, and N. K. Jha, “Sharks: Smart hacking approaches for risk scanning in internet-of-things and cyber-physical systems based on machine learning,” _IEEE Transactions on Emerging Topics in Computing_ , 2021.
  * [8] J. Brown, T. Saha, and N. K. Jha, “Gravitas: Graphical reticulated attack vectors for internet-of-things aggregate security,” _IEEE Transactions on Emerging Topics in Computing_ , 2021.
  * [9] D. P. Kingma and M. Welling, “Auto-encoding variational bayes,” _arXiv preprint arXiv:1312.6114_ , 2013.
  * [10] X. Zhang, S. Ji, and T. Wang, “On distinguishability criteria for estimating generative models,” _arXiv preprint arXiv:1412.6515_ , 2014.
  * [11] T. Salimans, I. Goodfellow, W. Zaremba, V. Cheung, A. Radford, and X. Chen, “Improved techniques for training gans,” in _NIPS_ , 2016.
  * [12] M. Arjovsky and L. Bottou, “Towards principled methods for training generative adversarial networks,” _arXiv preprint arXiv:1701.04862_ , 2017.
  * [13] T. Che, Y. Li, A. P. Jacob, Y. Bengio, and W. Li, “Mode regularized generative adversarial networks,” _arXiv preprint arXiv:1612.02136_ , 2016.
  * [14] X. Chen, Y. Duan, R. Houthooft, J. Schulman, I. Sutskever, and P. Abbeel, “Infogan: Interpretable representation learning by information maximizing generative adversarial nets,” in _NIPS_ , 2016.
  * [15] L. Metz, B. Poole, D. Pfau, and J. Sohl-Dickstein, “Unrolled generative adversarial networks,” _arXiv preprint arXiv:1611.02163_ , 2016.
  * [16] I. P. Durugkar, I. Gemp, and S. Mahadevan, “Generative multi-adversarial networks,” _CoRR_ , vol. abs/1611.01673, 2016.
  * [17] F. Huszar, “How (not) to train your generative model: Scheduled sampling, likelihood, adversary?” _CoRR_ , vol. abs/1511.05101, 2015.
  * [18] L. Theis, A. van den Oord, and M. Bethge, “A note on the evaluation of generative models,” _CoRR_ , vol. abs/1511.01844, 2015.
  * [19] I. J. Goodfellow, “Nips 2016 tutorial: Generative adversarial networks,” 2016\.
  * [20] A. Ghosh, V. Kulharia, V. P. Namboodiri, P. H. S. Torr, and P. K. Dokania, “Multi-agent diverse generative adversarial networks,” _CoRR_ , vol. abs/1704.02906, 2017.
  * [21] A. Ghosh, V. Kulharia, and V. P. Namboodiri, “Message passing multi-agent gans,” _CoRR_ , vol. abs/1612.01294, 2016. [Online]. Available: http://arxiv.org/abs/1612.01294
  * [22] M. Liu and O. Tuzel, “Coupled generative adversarial networks,” in _NIPS_ , 2016.
  * [23] Q. Hoang, T. D. Nguyen, T. Le, and D. Q. Phung, “Multi-generator generative adversarial nets,” _CoRR_ , vol. abs/1708.02556, 2017. [Online]. Available: http://arxiv.org/abs/1708.02556
  * [24] A. Odena, C. Olah, and J. Shlens, “Conditional image synthesis with auxiliary classifier gans,” in _ICML_ , 2017.
  * [25] D. J. Im, H. Ma, C. D. Kim, and G. W. Taylor, “Generative adversarial parallelization,” _CoRR_ , vol. abs/1612.04021, 2016. [Online]. Available: http://arxiv.org/abs/1612.04021
  * [26] A. Neelakantan, L. Vilnis, Q. V. Le, I. Sutskever, L. Kaiser, K. Kurach, and J. Martens, “Adding gradient noise improves learning for very deep networks,” _arXiv preprint arXiv:1511.06807_ , 2015.