Feature Ranking of Active Region Source Properties in Solar Flare Forecasting and the Uncompromised Stochasticity of Flare Occurrence

Solar flares are the most explosive events in the heliosphere, releasing in an abrupt way up to 10³³ erg of energy in a time interval typically ranging between 10 and 1000 s (Benz 2017). This energy is previously stored in specific magnetic configurations and, when magnetic reconnection occurs (Priest & Forbes 2002), it is transformed into mass acceleration, heating, and electromagnetic radiation at all wavelengths. It is also established that flares are a major space weather agent in the heliosphere (Schwenn 2006), while, as secondary effects through their correlation with coronal mass ejections, they induce geospace and ionospheric disturbances, malfunctions, and impairments on technologies in the geosphere, such as flight navigation, satellite communication, and power grid distribution.

Solar flare forecasting is a prominent discipline (Gallagher et al. 2002; Georgoulis & Rust 2007; Li et al. 2007, 2008; Schrijver 2007; Barnes & Leka 2008; Wang et al. 2008; Colak & Qahwaji 2009; Yu et al. 2009; Ahmed et al. 2013; Bobra & Couvidat 2015; Barnes et al. 2016; McCloskey et al. 2017; Murray et al. 2017; Sadykov & Kosovichev 2017; Benvenuto et al. 2018; Huang et al. 2018; Massone et al. 2018; Nishizuka et al. 2018; Park et al. 2018) within the recent field of space weather forecasting that relies on the availability of two ingredients; one observational and one computational. First, it is well-established that solar active regions (ARs) exclusively host major flares and therefore flare prediction needs experimental data on AR properties, associated to the photospheric and coronal magnetic field; however, coronal information has only recently started being used in the form of EUV images given as input to a deep learning network by Nishizuka et al. (2018). Second, this information on AR magnetic properties can be processed for prediction purposes by means of a computational method for data analysis; machine learning has recently offered strong candidates for such methods.

Since February 2010, the Helioseismic and Magnetic Imager on board the Solar Dynamics Observatory (SDO/HMI; Scherrer et al. 2012) is providing both line-of-sight and vector magnetograms of the full solar disk at a (vector magnetogram) cadence of 12 minutes. SDO/HMI magnetograms can be used for solar flare prediction according to two different approaches. First, HMI magnetograms are utilized to calculate a variety of properties either from the line-of-sight component only, from the radial component only, or from all three vector components. Various single-valued quantities, hereafter referred to as features, can be calculated from these property images through a variety of techniques (e.g., thresholding, feature recognition, etc.), such that calculation of one physical property may provide multiple features as inputs to machine learning (i.e., image maximum, total, and moments). Of course, additional features that are not derived from property images may also contribute to the input data set. Second, from a deep learning perspective, HMI images can be given as input to Convolutional Neural Networks that automatically perform a probabilistic forecasting. This present paper follows the first approach, and this is for several reasons. First, we had at our disposal the property extraction power provided by the algorithms developed within the Horizon 2020 FLARECAST project (http://flarecast.eu), which generated data sets of almost 200 features determined from properties extracted from photospheric SDO/HMI vector magnetograms. This database of features probably represents the highest data dimensionality currently available for flare forecasting purposes. Second, one of the objectives of our research was to determine to what extent AR properties are redundant when forecasting flares, and a straightforward way to do this is by ranking the extracted features according to their predictive capability. Finally, so far most publications in flare prediction utilize feature-based machine-learning methods, so another objective of this paper is to investigate how data preparation (and, specifically, the preparation of the training set in the case of supervised algorithms) impacts the prediction scores. Specifically, the analysis performed in this paper relies on two supervised machine-learning algorithms that combine prediction with feature ranking, namely hybrid LASSO (Benvenuto et al. 2018) and Random Forest (RF; Breiman 2001). However, two other methods of this kind, namely logit (Wu et al. 2009) and a support vector machine for classification (Cortes & Vapnik 1995), have been also applied to the same data sets for verification purposes, with coherent results.

The content of the paper is as follows. Section 2 overviews the data analysis procedure, describing in detail the features used for prediction, the data preparation process, and key aspects of the machine-learning methods adopted. Section 3 contains the results of the analysis, while Section 4 discusses these results. Our conclusions are offered in Section 5.

2.1. Data and Features

2.2. Data Preparation

2.3. Prediction and Feature Ranking Methods

The effectiveness of the prediction was assessed by skill scores computed on the previously unseen test sets. Following suggestions in Bloomfield et al. (2012), we chose to use the true skill statistic (TSS) and the Heidke skill score (HSS), assuming them as representative among skill scores existing in the literature. This said, although not shown here, we performed the analysis using the false alarm ratio, probability of detection, and accuracy metrics, obtaining similar results in terms of the relative forecasting effectiveness of the two machine-learning methods. We point out that all these scores are computed by means of binary predictions applied to the test set. However, as noted in Section 2.3, LASSO and RF are regression methods providing as outcome real positive numbers that can be interpreted in a probabilistic sense and, in our approach, the transformation of these variables into dichotomous yes/no responses is accomplished by applying a fuzzy clustering technique against the regression outcomes. Other works typically apply an arbitrary probability threshold, P_th, of 0.5 to create dichotomous forecasts (Leka et al. 2019a, 2019b), although it should be noted that discriminating thresholds optimized on TSS should find P_th values close to the climatology (i.e., the average flare-day rate; Bloomfield et al. 2012; Barnes et al. 2016). Figure 1 shows a comparison between the thresholding performances of the clustering technique and the ones provided by the optimization of TSS and HSS and by the use of an ROC curve. Interestingly, the two hybrid regression/clustering approaches provide similar results, which are rather conservative and rather close to the ones achieved by optimizing the HSS, especially in the case of the prediction of C1+ flares. Furthermore, the ROC curve method relies on cutoff values computed by means of the Youden index (Youden 1950), which formally leads to the maximization of the TSS; in fact the figure shows that the two values are always very close and the small differences are just related to the different numerical way the thresholding-search schemes were implemented. For C1+ flares, our hybrid approach results in P_th ≈ 0.4 for both HLA and RF, meaning that our C1+ TSS values are (probably) more comparable to those whose probabilities are converted to dichotomous forecasts using P_th = 0.5 (as our P_th lies closer to 0.5 than the average C1+ flare-day rate of ∼0.26). The situation is more complex for M1+ flares, however, as the average P_th found by the fuzzy-clustered HLA method is almost equivalent to that optimized on TSS (i.e., approaching the average flare-day rate of 0.04) while for the fuzzy-clustered RF method it instead occupies greater values that lie between the TSS and HSS optimized cases.

Zoom In Zoom Out Reset image size

The averages and standard deviations of the TSS and HSS values over 100 random realizations of the training/test sets for both prediction methods are shown in Table 1. In the case of HSS, the reliability of average values may be challenged by inconsistent flare/no-flare imbalance ratio across the 100 realizations. However, we have a posteriori checked the sample statistics of the 100 random realizations: average flare/no-flare imbalance ratios across the 100 test sets are ∼0.34 for C1+ events and ∼0.04 for M1+ events, with relative standard deviations that are <16% and <27% of these values, respectively (reflecting the largest relative standard deviations for the four separate UT issuing times considered here).

Focusing then on the feature ranking process, the boxplots in Figures 2 and 3 show the top 10 features ordered by their mean RFE ranking, obtained by HLA and RF over the 100 random realizations for each of the four forecast issuing times. Specifically, Figure 2 refers to the prediction of GOES C1+ flares, while Figure 3 refers to the prediction of GOES M1+ flares.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From these results it becomes possible to assess the impact of feature selection on the prediction performance, by computing specific skill scores and statistics in a cumulative way. The panels in Figures 4 and 5 plot the TSS values obtained by HLA and RF in the case of one specific data set realization, while adding one feature at a time, starting from the feature with the highest ranking, down to the feature with the tenth highest ranking. A given feature has the same color throughout each set of plots, for all issuing times.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to have a clearer picture of the features that repeatedly show the highest predictive impact, the histograms in Figure 6 compute the number of times over the four issuing times that each feature appears in the top-10 ranking of training set averages. These plots only present features that reach the top-10 ranking at least twice out of the four issuing times for a given machine-learning method and flare class. For the C1+ flare prediction (Figure 6; top row) one sees, for example, that the past flare history (flare_index_past) and Schrijver's (Schrijver 2007) R-value (r_value_br_logr) consistently appear for both HLA and RF. This is not the case for the prediction of M1+ flares (Figure 6; bottom row). It should be noted that the importance of the specific features may only be due to the machine-learning method used; it is their consistency of appearance, however, that is notable.

Zoom In Zoom Out Reset image size

We first notice that the maximum values of HSS and TSS achieved in Table 1 are distinctly different from one, indicating far from perfect performance. Interestingly, these scores are almost systematically smaller than the ones recently achieved by methods illustrated in Bobra & Couvidat (2015) and Florios et al. (2018) that use input data with a significantly smaller dimensionality. The methods described in those papers are all supervised, utilize features extracted from HMI data, and perform predictions in a 24 h window. However, the way data preparation is performed and, in particular, how the training set is constructed is significantly different than what is done in the present paper. In particular:

• 1.  
  The test sets utilized in this work to assess the performances of HLA and RF do not contain feature vectors belonging to ARs with feature vectors contained in the training sets. Instead, the training sets utilized in Bobra & Couvidat (2015), and Florios et al. (2018) combined feature vectors belonging to the same ARs in the two sets.
• 2.  
  We constructed four separate training/testing sets, each corresponding to a specific UT forecast issuing time on all of the days considered. Our results show reasonably consistent forecast performance across these four issuing-time sets. However, the main benefit to this approach is in the interpretation of the feature selection results. Identifying key forecasting features through their appearance in all (or most) of the top-10 feature ranking lists across these four issuing-time sets increases their robustness through temporal consistency.
• 3.  
  The training set utilized in Florios et al. (2018) is populated with approximately the same number of vectors as the test set, while in our approach (and in the one followed in Bobra & Couvidat 2015) the machine-learning methods are trained with training sets two times more populated than the test sets, which is more realistic with respect to typical experimental settings.
• 4.  
  Our prediction methods are optimized using a fuzzy clustering technique, while in the other three cases the input parameters are fixed in such a way to optimize a specific skill score (namely, the TSS).

In order to assess the impact of these differences against the methods' performance, we trained HLA and RF using all 14,931 point-in-time feature vectors distributing them between training and test sets as was done in Bobra & Couvidat (2015) and Florios et al. (2018). Specifically, we generated the training and test sets focusing on feature vectors instead of on ARs, i.e., we randomly extracted the feature vectors from the database at disposal without imposing any constraint that forbids feature vectors of the same AR to populate both the training and the test set (the two sets are populated as in Bobra & Couvidat 2015, following a 2:1 proportion). Furthermore, we did not care for time consistency and so we mixed up feature vectors belonging to different issuing times. Finally, the prediction methods are optimized in such a way to maximize the TSS. Table 2 shows that TSS increases significantly in the prediction of C1+ and M1+ flares for both HLA and RF and also HSS produced by RF becomes larger, although less significantly. These scores are now more in line with the ones obtained in the other two papers, at the same time showing smaller standard deviations. This leads to the conclusion that, not surprisingly, also in flare prediction the biases introduced in the process of training set generation and the way the algorithms are optimized strongly influence the performance of supervised methods.

As far as feature ranking is concerned, it is evident from Figures 2 and 3 that first, features with the best ranking have the smallest standard deviations, so their impact on prediction is consistently high, regardless of the splitting of the data used for training the machine-learning algorithms. In the case of prediction of GOES C1+ flares, colors largely repeat in all four panels, telling us that features with highest predictive power are common to all issuing times considered. This behavior is not as robust in the case of prediction of GOES M1+ flares, but we are confident that this is a consequence of the lower occurrence rate of M1+ flares and the resulting variation in the flare/no-flare imbalance ratio of the random training sets. A consistent imbalance ratio is more or less guaranteed for C1+ flares, whose comprehensive statistics over solar cycle 24 ensure a well-balanced training process.

Figures 4 and 5 show that only a small number of features (up to 10) over the scores of features proposed and/or applied for flare prediction, are sufficient to achieve maximum performance of a given machine-learning method. Notice from these figures that the highest-ranking feature alone (feature 1) suffices to give TSS and HSS values that are at least half of the maxima achieved. Up to the fourth feature, the values of TSS and HSS saturate already, indicating that adding more features will not improve (and may, in fact, be detrimental to) prediction performance. Also, provided that flare statistics are sufficient to deal with the flare/no-flare imbalance ratio in the random selection of training and test sets, these few best-performing features are consistent for a given prediction method. In their study, Bobra & Couvidat (2015) found that the four most significant features in their analysis were the total unsigned current helicity, total magnitude of the Lorentz force, total photospheric magnetic free energy density, and unsigned vertical current. In our study, Figure 6 shows that features associated with the unsigned vertical current (i.e., total, maximum, standard deviation) are among the most temporally consistent of our best-performing features, particularly for C1+ flares (its standard deviation is in the top-ten features of all four UT issuing times for Hybrid LASSO, while its total and maximum are in the top-10 of three of the four UT issuing times for RF) and less so for M1+ flares (its standard deviation is in the top-10 for two of the four UT issuing times). However, the same figure shows that these predictors may well change from method to method, which hampers efforts to understand physically why some features work better than others. Best performers appear to change also for the prediction of different flare classes, which is a very interesting finding that undoubtedly warrants additional investigation in the future.

This study employs the highest-dimension data set of prediction features to date in regards to solar flare forecasting, while it shows TSS values similar to the better performing region-by-region forecasting systems in the literature. This point taken, the actual HSS and TSS values are not identical to—and may even be somewhat lower than—the respective values reported in Bobra & Couvidat (2015) and in Florios et al. (2018), with the latter study using RF applied to FLARECAST data, as well. The reason is that training and testing of machine-learning methods in the present study were not only performed on nonoverlapping data as in previous studies, but even the solar ARs selected for training and testing were different. This conclusion seems in line with the considerations contained in Camporeale (2019), whose careful assessment of the causality issue identifies this as one of the crucial aspects impacting the forecasting performances.

The rationale for using hybrid LASSO and RF in this work is these methods' ability to perform feature ranking via RFE, among other methods (i.e., Fisher's score, Gini index, etc.). However, there are 26 machine-learning methods implemented in FLARECAST. Their definitive evaluation is in progress, so the values of pertinent skill scores may well increase in future studies utilizing FLARECAST data, in the search for finding the optimal machine-learning method(s) for the near-realtime FLARECAST forecasting service. We also understand that a meaningful methodological next step would be to introduce deep learning methods in the pipeline. However, interestingly, the use of these more modern approaches in flare forecasting does not necessarily imply significantly higher skill scores (see, e.g., (Nishizuka et al. 2018), where TSS and HSS for the prediction of M1+ flares are reported as 0.80 and 0.26, respectively).

What will, most likely, not change in the foreseeable future are the following two core conclusions of this work.

First, the current range of properties that have been extracted from the HMI magnetograms show significant redundancy and no more than 10 features contained in these properties are sufficient to allow machine-learning methods to achieve maximum performance.

Second, and perhaps foremost, the maximum values of HSS and TSS achieved are distinctly different from one, indicating far from perfect performance. In physical terms, even using the largest flare prediction data volume assembled to date, we have not managed to substantially surpass the performance of a random-chance forecast (as shown by HSS) or to substantially increase the probability of detection (0.57–0.65 for C1+) despite an encouragingly low probability of false detection (∼0.10). The latter two compete against each other to result in TSS. This is equivalent to saying that flare prediction remains probabilistic, rather than binary yes/no with a perfect performance. The core reasons for this may be multiple: first, we only rely on photospheric magnetic field data, but flares occur above the line-tied photosphere in the low solar corona. Second, flares may well be intrinsically stochastic phenomena, as adopted in a long-standing working hypothesis (Rosner & Vaiana 1978), shown conclusively by the flares' time-dependent Poisson waiting times (Crosby et al. 1998; Wheatland & Litvinenko 2002) and interpreted physically via the concept of self-organized criticality (Lu & Hamilton 1991; Lu et al. 1993; Vlahos et al. 1995)—see also Aschwanden et al. (2016) for a comprehensive review.

This study was enabled by the European Union's Horizon 2020 Research and Innovation Action Flare Likelihood And Region Eruption foreCASTing (FLARECAST) project, under grant agreement No. 640216. All authors warmly acknowledge the support of the FLARECAST project.

Here we provide details of the FLARECAST feature labels used in this work, with short descriptions and references to their original definition/implementation (or, e.g., detection methods used in their calculation). Features are grouped in the following manner: Table 3 contains those features derived from B_los only and B_radial only (Georgoulis 2005, 2012; Georgoulis & Rust 2007; Schrijver 2007; Falconer et al. 2008; Hewett et al. 2008; Ahmed et al. 2010; Mason & Hoeksema 2010; Georgoulis et al. 2012; Zuccarello et al. 2014; Guerra et al. 2015; Kontogiannis et al. 2018); Table 4 contains those features requiring all three vectormagnetic field components (Kusano et al. 2002; Schuck 2008); Table 5 contains only those features related to the total and mean quantities provided as the SHARP keywords of Bobra et al. (2014).

Table 3.  FLARECAST B_los- and B_r-derived Feature List with Short Descriptions

Feature Label	Description
alpha_exp_fft_blos/alpha, alpha_exp_fft_br/alpha	Fourier power spectral index
alpha_exp_cwt_blos/alpha, alpha_exp_cwt_br/alpha	Continuous wavelet transform power spectral index
beff_blos/beff, beff_br/beff	Effective connected magnetic field strength Beff
decay_index_blos/max_l_over_hmin	Max. ratio of MPIL length to min. height of critical decay index
decay_index_br/max_l_over_hmin	l/h(ncr)min
decay_index_blos/tot_l_over_hmin, decay_index_br/tot_l_over_hmin	Total of all separate MPIL ratios of l/h(ncr)min
decay_index_blos/l_over_minhmin, decay_index_br/l_over_minhmin	Ratio of MPIL l/h(ncr)min (for MPIL having lowest h(ncr)min)
decay_index_blos/maxl_over_hmin, decay_index_br/maxl_over_hmin	Ratio of MPIL l/h(ncr)min (for longest MPIL)
flare_past	Binary flag for occurrence of ≥1 flare in previous 24 hr
flare_index_past	Accumulated GOES flare peak magnitudes in previous 24 hr
frdim_blos/frdim, frdim_br/frdim	Fractal dimension
gs_slf/g_s	Sum of the horizontal magnetic gradient
gs_slf/slf	Separation distance lead. and follow. polarity subgroups
ising_energy_blos/ising_energy, ising_energy_br/ising_energy	Ising energy (calculated pixel-by-pixel)
ising_energy_part_blos/ising_energy_part, ising_energy_part_br/ising_energy_part	Ising energy (calculated using Beff flux partitions)
lat_hg	Heliographic latitude of SHARP centroid
lon_hg	Heliographic longitude of SHARP centroid
mf_spectrum_blos/dq, mf_spectrum_br/dq	Multifractal generalized correlation dimension spectrum
mpil_blos/max_length, mpil_br/max_length	Maximum length of a single MPIL
mpil_blos/tot_length, mpil_br/tot_length	Total length of all MPILs
mpil_blos/tot_usflux, mpil_br/tot_usflux	Total unsigned flux around all MPILs
r_value_blos_logr, r_value_br_logr	Schrijver's R (log10 form)
sfunction_blos/zq, sfunction_br/zq	Multifractal structure function inertial range index
wlsg_blos/value_int, wlsg_br/value_int	Falconer's (${}^{{\rm{L}}}{\mathrm{WL}}_{\mathrm{SG}}$

Download table as:  ASCIITypeset image