3.3 Nonlinear methods

“...every complicated question has a simple answer which is wrong. Analyzing atime series with a nonlinear approach is definitely a complicated problem. Simpleanswers have been repeatedly offered in the literature, quoting numerical valuesfor attractor dimensions for any conceivable system.” (Hegger
et al., 1999)

The nonlinearities in the dynamo equations readily give rise to chaotic behaviour of the
solutions. The long term behaviour of solar activity, with phenomena like grand minima and grand
maxima, is also suggestive of a chaotic system. While chaotic systems are inherently unpredictable
on long enough time scales, their deterministic nature does admit forecast within a limited
range. It is therefore natural to explore this possibility from the point of view of solar cycle
prediction.

3.3.1 Attractor analysis and phase space reconstruction: the pros ...

Assuming that the previous values of the sunspot number do in some way determine the current
expected value, our problem becomes restricted to an -dimensional phase space, the dimensions being
the current value and the previous values. With a time series of length , we have
points fixed in the phase space, consecutive points being connected by a line. This
phase space trajectory is a sampling of the attractor of the physical system underlying the
solar cycle (with some random noise added to it). The attractor represents a mapping in phase
space which maps each point into the one the system occupies in the next time step: if this
mapping is known to a good approximation, it can be used to extend the trajectory towards the
future.

For the mapping to be known, needs to be high enough to avoid self-crossings in the phase space
trajectory (otherwise the mapping is not unique) but low enough so that the trajectory still yields a good
sampling of the attractor. The lowest integer dimension satisfying these conditions is the embeddingdimension of the attractor (which may have a fractal dimension itself).

Once the attractor has been identified, its mathematical description may be done in two
ways.

(1) Parametric fitting of the attractor mapping in phase space. The simplest method is the piecewise
linear fit suggested by Farmer and Sidorowich (1987) and applied in several solar prediction attempts, e.g.,
Kurths and Ruzmaikin (1990). Using a method belonging to this group, for cycle 24 Kilcik et al. (2009)
predict a peak amplitude of 87 to be reached in late 2012. Alternatively, a global nonlinear fit can also be
used: this is the method applied by Serre and Nesme-Ribes (2000) as the first step in their global flow
reconstruction (GFR) approach.

(2) Nonparametric fitting. The simplest nonparametric fit is to find the closest known attractor point to
ours (in the -dimensional subspace excluding the last value) and then using this for a prediction,
as done by Jensen (1993). (This resulted in so large random forecast errors that it is practically unsuitable
for prediction.) Neural networks, discussed in more detail in Section 3.3.4 below, are a much more
sophisticated nonparametric fitting device.

(3) Indirectly, one may try to find a set of differential equations describing a system that gives rise to an
attractor with properties similar to the observed. In this case there is no guarantee that the derived
equations will be unique, as an alternative, completely different set may also give rise to a very
similar attractor. This arbitrariness of the choice is not necessarily a problem from the point
of view of prediction as it is only the mapping (the attractor structure) that matters. Such
phase space reconstruction by a set of governing equations was performed, e.g., by Serre and
Nesme-Ribes (2000) or Aguirre et al. (2008); for cycle 24 the latter authors predict a peak
amplitude of 65 ± 16. On the other hand, instead of putting up with any arbitrary set of
equations correctly reproducing the phase space, one might make an effort to find a set with a
structure reasonably similar to the dynamo equations so they can be given a meaningful physical
interpretation. Methods following this latter approach will be discussed in Sections 4.4 and
4.5.

3.3.2 ... the cons ...

Finding the embedding dimension and the attractor structure is not a trivial task, as shown by the widely
diverging results different researchers arrived at. One way to find the correct embedding dimension is the
false nearest neighbours method (Kennel et al., 1992), essentially designed to identify self-crossing in the
phase space trajectory, in which case the dimension needs to be increased. But self-crossings are to
some extent inevitable, due to the stochastic component superimposed on the deterministic skeleton of the
system.

As a result, the determination of the minimal necessary embedding dimension is usually done indirectly.
One indirect method fairly popular in the solar community is the approach proposed by Sugihara and
May (1990) where the correct dimension is basically figured out on the basis of how successfully the model,
fit to the first part of the data set, can “predict” the second part (using a piecewise linear
mapping).

Another widely used approach, due to Grassberger and Procaccia (1983), starts by determining the
correlation dimension of the attractor, by simply counting how the number of neighbours in an embedding
space of dimension increases with the distance from a point. If the attractor is a lower
dimensional manifold in the embedding space and it is sufficiently densely sampled by our data then the
logarithmic steepness of this function should be constant over a considerable stretch of the curve: this is
the correlation dimension . Now, we can increase gradually and see at what value saturates:
that value determines the attractor dimension, while the value of where saturation is reached yields
the embedding dimension.

The first nonlinear time series studies of solar activity indicators suggested a time series spacing of
2 – 5 years, an attractor dimension 2 – 3 and an embedding dimension of 3 – 4 (Kurths and
Ruzmaikin, 1990; Gizzatullina et al., 1990). Other researchers, however, were unable to confirm these
results, either reporting very different values or not finding any evidence for a low dimensional attractor at
all (Calvo et al., 1995; Price et al., 1992; Carbonell et al., 1994; Kilcik et al., 2009; Hanslmeier and
Brajša, 2010). In particular, I would like to call attention to the paper by Jensen (1993), which, according
to ADS and WoS, has received a grand total of zero citations (!) up to 2010, yet it displays an exemplary
no-nonsense approach to the problem of sunspot number prediction by nonlinear time series methods.
Unlike so many other researchers, the author of that paper does not fool himself into believing to see a
linear segment on the logarithmic correlation integral curve (his Figure 4); instead, he demonstrates on a
simple example that the actual curve can be perfectly well reproduced by a simple stochastic
process.

These contradictory results obviously do not imply that the mechanism generating solar activity is not
chaotic. For a reliable determination a long time series is desirable to ensure a sufficiently large number of
neighbours in a phase space volume small enough compared to the global scale of the attractor. Solar data
sets (even the cosmogenic radionuclide proxies extending over millennia but providing only a decadal
sampling) are typically too short and sparse for this. In addition, clearly distinguishing between the phase
space fingerprints of chaotic and stochastic processes is an unsolved problem of nonlinear dynamics
which is not unique to solar physics. A number of methods have been suggested to identify
chaos unambiguously in a time series but none of them has been generally accepted and this
topic is currently a subject of ongoing research – see, e.g., the work of Freitas et al. (2009)
which demonstrates that the method of “noise titration”, somewhat akin to the Sugihara–May
algorithm, is uncapable of distinguishing superimposed coloured noise from intrinsically chaotic
systems.

3.3.3 ... and the upshot

Starting from the 1980s many researchers jumped on the chaos bandwagon, applying nonlinear time series
methods designed for the study of chaotic systems to a wide variety of empirical data, including solar
activity parameters. From the 1990s, however, especially after the publication of the influential book by
Kantz and Schreiber (1997), it was increasingly realized that the applicability of these nonlinear algorithms
does not in itself prove the predominantly chaotic nature of the system considered. In particular, stochastic
noise superposed on a simple, regular, deterministic skeleton can also give rise to phase space characteristics
that are hard to tell from low dimensional chaos, especially if strong smoothing is applied to
the data. As a result, the pendulum has swung in the opposite direction and currently the
prevailing view is that there is no clear cut evidence for chaos in solar activity data (Panchev and
Tsekov, 2007).

One might take the position that any forecast based on attractor analysis is only as good as the
underlying assumption of a chaotic system is: if that assumption is unverifiable from the data, prediction
attempts are pointless. This, however, is probably a too hasty judgment. As we will see, the
potentially most useful product of phase space reconstruction attempts is the inferences they allow
regarding the nature of the underlying physical system (chaotic or not), even offering a chance to
constrain the form of the dynamo equations relevant for the Sun. As discussed in the next
section, such truncated models may be used for forecast directly, or alternatively, the insight they
yield into the mechanisms of the dynamo may be used to construct more sophisticated dynamo
models.

3.3.4 Neural networks

Neural networks are algorithms built up from a large number of small interconnected units (“neurons” or
“threshold logic units”), each of which is only capable of performing a simple nonlinear operation on an
input signal, essentially described by a step function or its generalized (rounded) version, a sigmoid
function. To identify the optimal values of thresholds and weights parameterizing the sigmoid functions of
each neuron, an algorithm called “back propagation rule” is employed which minimizes (with or
without human guidance) the error between the predicted and observed values in a process called
“training” of the network. Once the network has been correctly trained, it is capable of further
predictions.

The point is that any arbitrary multidimensional nonlinear mapping may be approximated by a
combination of stepfunctions to a good degree – so, as mentioned in Section 3.3.1 above, the neural network
can be used to find the nonlinear mapping corresponding to the attractor of the given time
series.

More detailed introductions to the method are given by Blais and Mertz (2001) and by Calvo
et al. (1995); the latter authors were also the first to apply a neural network for sunspot number prediction.
Unfortunately, despite their claim of being able to “predict” (i.e., postdict) some earlier cycles
correctly, their prediction for cycle 23 was off by a wide margin (predicted peak amplitude of
166 as opposed to 121 observed). One of the neural network forecasts for cycle 24 (Maris and
Oncica, 2006) was equally far off, predicting a maximum of 145 as early as December 2009, while
another one (Uwamahoro et al., 2009) yields a more conservative value of 117.5 ± 8.5 for
2012.