### 3.3 Nonlinear methods

“...every complicated question has a simple answer which is wrong. Analyzing a time series with a nonlinear approach is definitely a complicated problem. Simple answers have been repeatedly offered in the literature, quoting numerical values for 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 embedding dimension 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, 1990Gizzatullina 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., 1995Price et al., 1992Carbonell et al., 1994Kilcik et al., 2009Hanslmeier 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.