For response variables, we take advantage of the recorded log intensities of flare events in the GOES data set (Garcia, 1994) ranging from 2 June 2010 to 29 December 2018. The flare events are recorded at their peak time (time at the highest flare intensity). Although the theoretical distribution of the flare events should be a power law distribution, the observed distribution is different from the theoretical distribution because flares in lower-energy levels are lost in the background and go undetected (Jiao et al., 2020). In this paper, we focus on the observed information. By scientific convention, solar flares belong to category B if their log intensity (log₁₀) is within (−∞, −6), category C if [−6, −5), category M if [−5, −4), and category X if (−4, ∞). Figure 1D shows that the M/X flare events are far fewer than B/C. The data imbalance issue is addressed in Section 3.1.

For covariates/features, we consider SHARP data (Schou et al., 2012) from 860 HARPs during the same time period (2 June 2010 to 29 December 2018). Approximately 7,000 HARPs are observed, many occurring without flares. From these, to maintain the quality of the data, we down-select the HARPs to a group of 860 based on the criteria that 1) the longitude of the HARP should be within the range of ±68° from the central meridian of the Sun to avoid projection effects (Bobra and Couvidat, 2015; Chen et al., 2019) and 2) the missing SHARP parameters should be fewer than 5% of all in the HARP to make sure that the missing data are not significantly large to cause any bias in model training.

For this type of data, an important practical goal of any model is to forecast the future flare intensity, given an observed value of SHARP parameters. As such, for each point in our dataset, we match the corresponding SHARP covariates with the GOES flare list (response variable) at the time point that is equal to the peak time subtracted by Δt, where Δt is the prediction time window and can take values in {6, 12, 24, 36, 48} h.

All covariates in the raw data are shown in Table 2. These are the same features used by Chen et al. (2019) and Bobra and Couvidat (2015). However, we only use a subset because our correlation analysis showed that some features are extremely highly correlated. Specifically, TOTUSJH/TOTUSJZ (0.9959), SIZE_ACR/NACR (0.9999), SHRGT45/MEANSHR (0.9969), and SIZE/NPIX (0.9999) are extremely highly correlated. For each of these features, we keep one feature and leave out the other. After removal, we are left with 16 covariates: TOTUSJH, SAVNCPP, USFLUX, ABSNJZH, TOTPOT, SIZE_ACR, MEANPOT, SIZE, MEANJZH, SHRGT45, MEANJZD, MEANALP, MEANGBT, MEANGAM, MEANGBZ, and MEANGBH. Their correlation matrix is displayed in Figure 1B.

A key aspect of our models is modeling the flare intensities by ARs. Consequently, we need to ensure similar sets of ARs in both training and testing. To facilitate this, we randomly split data into a training set and testing set, the latter containing 20% of the total across all the ARs with two or more events. For those with only one event, we flip a biased coin with probability 0.2, assigning to the testing set if the result was true and, otherwise, to the training set. This random split scheme guarantees a fair representation of each AR in both the training set and testing set. Each feature in the training set was then standardized to have a zero mean and a standard deviation of one. We used these parameters to normalize the testing set. Additionally, we further randomly divided the training set by ARs into a sub-train set and a validation set, with the latter equal to 20% of the original set. We used this validation set for the model selection process discussed in Section 4.

Weighted likelihood and weighted maximum likelihood estimation (MLE) have been used in the literature for robust estimations, especially when outliers exist in the data (Carroll and Pederson, 1993; Field and Smith, 1994; Markatou et al., 1998). The idea is to down-weight the outlier data points so that they do not deteriorate the performance of the model too much. A similar principle can be applied here to handle the imbalance problem, i.e., down-weighting the “majority” data points and/or up-weighting the “minority” points.

To illustrate, we consider a simple example of a standard linear regression setting. Assume that the response y _i given covariate X _i follows the normal distribution y i | X i ∼ N ( X i T β , 1 ) , i = 1 , … , n , where β is the linear regression coefficient. The least squares estimate for β is given by β ̂ = arg max β ∑ i = 1 n − ( y i − X i T β ) 2 . With highly imbalanced data like the solar flare data given in Section 2, standard linear regression would not work well, and weighted linear regression may be used instead. In that case, we need to calculate the weighted estimator β ̃ = arg max β ∑ i = 1 n − w i ⋅ ( y i − X i T β ) 2 for known weights w _i that may depend on the data. Moving beyond this simple example, the modeling we propose in the next sections is based on mixture models. Consider a simplified generative model as follows: y i | X i , z i = k ∼ N y i | X i T β k , σ k 2 , k = 1 , … , K , i = 1 , … , n , where K denotes the number of mixture components and z _i is a categorical variable with P (z _i = k) = π _k for all k; ∑ _k π _k = 1. Applying the weighted likelihood idea, we aim to find the optimizer arg max θ ∑ i = 1 n w i ⋅ log ⁡ p ( y i | X i , θ ) , where θ ≔ { β k , σ k 2 , π k } 1 K . Recall that the original log-likelihood is l ( θ ) ≔ ∑ i = 1 n ⁡ log ⁡ p ( y i | X i , θ ) . For mixture models, it is not straightforward to directly maximize l(θ). The expectation-maximization (EM) algorithm is needed (Dempster et al., 1977). Typically, the EM algorithm works with the logarithm of the complete data likelihood, which is defined as l ( π , β , σ 2 ) ≔ ∑ k = 1 K ∑ i = 1 n 1 ( z i = k ) ⋅ log π k + log N ( y i | X i T β k , σ k 2 ) . In the E-step, the EM algorithm computes τ i , k ≔ p ( z i = k | X i , y i ) ≔ π k N ( y i | X i T β k , σ k 2 ) ∑ j = 1 K π j N ( y i | X i T β j , σ j 2 ) , k = 1, … , K; then, it maximizes the expected complete log-likelihood, arg max β , σ 2 , π Q ( β , σ 2 , π ) ≔ arg max β , σ 2 , π ∑ i = 1 n ∑ k = 1 K τ i , k ⋅ log π k + log N ( y i | X i T β k , σ k 2 ) , in the M-step. That yields π ̂ k = ∑ i = 1 n τ i k n and μ ̂ k = ∑ i = 1 n τ i , k X i ∑ i = 1 n τ i , k and σ ̂ k 2 = ∑ i = 1 n τ i k ⋅ ( X i − μ ̂ k ) 2 ∑ i = 1 n τ i k ; k = 1, …, K.

To adapt the EM framework to find the MLE for the weighted likelihood, we note that under the EM algorithm, the log-likelihood is lower-bounded by the expected complete log-likelihood, i.e., l(θ) ≥ Q(θ), and by optimizing the lower bound Q(θ), the EM algorithm yields the (local) maximum of the likelihood l(θ) (Bishop, 2006). Under the weighted likelihood setting, we can also find a lower bound as follows: ∑ i = 1 n w i ⁡ log ⁡ p y i | X i , θ = ∑ i = 1 n w i ⁡ log ∑ k = 1 K p y i , z i = k | X i , θ = ∑ i = 1 n w i ⋅ log ∑ k = 1 K q z i = k ⋅ p y i , z i = k | X i , θ q z i = k ≥ ∑ i = 1 n w i ⋅ ∑ k = 1 K q z i = k log p y i , z i = k | X i , θ q z i = k .

The last inequality is due to Jensen’s inequality. Since the logarithm is a concave function, and ∑ k = 1 k q ( z i = k ) = 1 , we moved the log inside and left q outside. Here, we can follow the usual EM procedure to find the optimal q(z). The E-step sets q ( z i = k ) = P ( z i = k | X i , y i ) , and then, the M-step maximizes the subsequent expression. Therefore, we observe that the weighted log-likelihood is bounded below by the weighted expected complete log-likelihood. In other words, to find the weighted log-likelihood estimator with the EM framework, we can optimize the lower bound ∑ i = 1 n ∑ k = 1 K w i ⋅ τ i , k ⋅ log π k + log N ( y i | X i T β k , σ k 2 ) . Moreover, since w _i (⋅) is a known weight function, if we ensure that 0 < w _i ≤ C < ∞ for all i for some constant C free of parameters, then, by noting that log  π _k < 0, we have ∑ i = 1 n ∑ k = 1 K w i ⋅ τ i , k ⋅ log π k + log N ( y i | X i T β k , σ k 2 ) = ∑ i = 1 n ∑ k = 1 K τ i , k ⋅ w i ⋅ log π k + w i ⁡ log N ( y i | X i T β k , σ k 2 ) ≥ ∑ i = 1 n ∑ k = 1 K τ i , k ⋅ C ⋅ log π k + w i ⁡ log N ( y i | X i T β k , σ k 2 ) .

Because ∑ k = 1 K π k = 1 , it is easy to derive that argmax  β , σ 2 , π ∑ i = 1 n ∑ k = 1 K τ i , k ⋅ C ⋅ log π k + w i ⁡ log N y i | X i T β k , σ k 2 = argmax  β , σ 2 , π ∑ i = 1 n ∑ k = 1 K τ i , k ⋅ log π k + w i ⁡ log N y i | X i T β k , σ k 2 .

Therefore, the above expression is the lower bound of the weighted complete log-likelihood, and optimizing it, in turn, increases the data likelihood. The final expression is the lower bound to be optimized for our mixture models, which is proposed in the subsequent sections. The justification of the algorithm comes from the fact that the EM-based inference algorithm will converge to a local maximum of the weighted likelihood function (Wu, 1983).

In this section, let w _i be denoted by w (y _i ) to emphasize the fact that the weights only depend on y _i (s). By scientific convention, the log intensity is (in log₁₀) y _i ∈ (−∞, −6) for flare category B, y _i ∈ [−6, −5) for flare category C, y _i ∈ [−5, −4) for flare category M, and y _i ∈ (−4, ∞) for flare category X. We propose the following scheme for w (y _i ): w y i = 1 ∑ i = 1 n 1 y i ≥ − 4 / n  if  y i ≥ − 4 1 ∑ i = 1 n 1 − 5 ≤ y i < − 4 / n  if  − 5 ≤ y i < − 4 1 ∑ i = 1 n 1 − 6 ≤ y i < − 5 / n  if  − 6 ≤ y i < − 5 1 ∑ i = 1 n 1 y i < − 6 / n  if  y i < − 6 .

Note that as n → ∞, by the strong law of large numbers, w y → a . s . 1 ∫ − 4 ∞ p 0 y d y  if  y ≥ − 4 1 ∫ − 5 − 4 p 0 y d y  if  − 5 ≤ y < − 4 1 ∫ − 6 − 5 p 0 y d y  if  − 6 ≤ y < − 5 1 ∫ − ∞ − 6 p 0 y d y  if  y < − 6 , where p ₀ y) is the marginal distribution of y. The justification for choosing w(y) this way is as follows. As the number of data points goes to infinity, and it is assumed that (X, y) follows the “true” distribution p ₀ (X, y), by the law of large numbers, the weighted score function becomes 1 n ∑ i = 1 n w y i log ⁡ p y i | X i , θ → a . s . ∫ w y log ⁡ p y | x , θ p 0 x , y d x d y = ∫ w y ⋅ ∫ log ⁡ p y | x , θ p 0 x | y d x ⋅ p 0 y d y .

r (y, θ)≔∫ log  p (y|x, θ) p ₀ (x|y)dx and I _B ≔(−∞, −6), I _C ≔[−6, −5), I _M ≔[−5, −4), I _X ≔[−4, ∞) are defined. These intervals correspond to the scientific thresholds of flare categories B, C, M, and X, respectively. As the number of data points goes to infinity, the weighted log-likelihood function can now be written as ∫ w ( y ) log ⁡ p ( y | x , θ ) p 0 ( x , y ) d x d y = ∫ w ( y ) r ( y , θ ) p 0 ( y ) d y = ∫ I B w ( y ) r ( y , θ ) p 0 ( y ) d y + ∫ I C w ( y ) r ( y , θ ) p 0 ( y ) d y + ∫ I M w ( y ) r ( y , θ ) p 0 ( y ) d y + ∫ I X w ( y ) r ( y , θ ) p 0 ( y ) d y .

As mentioned previously, the number of data points of M and X are fewer than that of B and C (Figure 1D), i.e., p ₀ y) places negligible masses on I _M and I _X compared to I _B and I _C . As a consequence, the four terms in the last expression are of different scales. As such, setting w(y) to the above-proposed choice effectively “normalizes” and places these four components on the same scale and balances them out. Explicitly, ∫ w ( y ) log ⁡ p ( y | x , θ ) p 0 ( x , y ) d x d y = ∫ I B r ( y , β ) p 0 ( y ) d y ∫ I B p 0 ( y ) d y + ∫ I C r ( y , β ) p 0 ( y ) d y ∫ I C p 0 ( y ) d y + ∫ I M r ( y , β ) p 0 ( y ) d y ∫ I M p 0 ( y ) d y + ∫ I X r ( y , β ) p 0 ( y ) d y ∫ I X p 0 ( y ) d y .

Finally, to tie this back to the last section, w(y) is then used as the weights in the weighted complete log-likelihood. For instance, in the example given in the previous section, the weighted complete log-likelihood is l ( π , β , σ 2 ) ≔ ∑ k = 1 K ∑ i = 1 n 1 ( z i = k ) ⋅ w ( y i ) ⋅ log π k + log N ( y i | X i T β k , σ k 2 ) .

As described in Section 2, X _i and y _i denote the SHARP parameter covariates and the corresponding log flare intensity response variable, respectively. We reiterate that the value of y _i is measured at time Δt h behind the values of X _i since we want to model y _i as the forecasted flare intensity at Δt h after observing SHARP values of X _i . Here, we use mixture modeling to characterize both the heterogeneity and shared patterns among active regions. Let K > 1 be the number of mixture components. We model the heterogeneity to be shared across ARs by the global parameters β _k and σ k 2 for k = 1, 2, … , K. For r = 1, … , R, each active region r is equipped with a discrete latent variable z ^r taking values in {1, 2, … , K}. If z ^r = k, then all the flare events { X i r , y i r } in the active region r follow the normal distribution y i r ∼ N ⋅ | β k T ⋅ X i r , σ k 2 . Latent variable z ^r is introduced to capture the “intrinsic” categories of flaring mechanisms from SDO/HMI data. The values of z ^r do not necessarily correspond to the scientific B/C/M/X categories but, rather, are the inferred “clusters” from the data. One important constraint under MM-R is that all the events under the same active region must have the same regression pattern parameterized by β z r , σ z r 2 .

Mathematically, the model is defined as follows. For AR r = 1, … , R, z r | π 1 : K ∼ Cat ⋅ | π y i r | z r = k ∼ N ⋅ | β k T X i r , σ k 2 , i = 1 , … , n r .

The model can also be represented as a probabilistic graphical model (Koller and Friedman, 2009), as shown in Figure 2. There are 2K + 1 “global” parameters { β k , σ k 2 } k = 1 K and π. Using the plate notation in this figure, there are R conditionally i.i.d. latent variables { z r } r = 1 R for each of the unique active regions. Under each active region r, there are n _r flare events { ( X i r , y i r ) } i = 1 n r .

Parameter estimation is achieved through an expectation-maximization algorithm (Dempster et al., 1977), with the lower bound as a weighted complete log-likelihood, as defined in section 3.2, and we optimize the weighted complete log-likelihood to remediate the data imbalance problem. Under this model, rather than the standard expected complete log-likelihood, our EM algorithm optimizes the weighted expected complete likelihood optimization problem, argmax π , σ 2 , β ∑ k = 1 K ∑ r = 1 R E z r = k | X r , y r log π k − ∑ i = 1 n r w i 2 log ⁡ σ k 2 + w i 2 σ k 2 ⋅ y i r − β k T X i r 2 .

The E-step computes τ k r = p z r = k | X r , y r = π k ⋅ ∏ i = 1 n r N y i r | β k T X i r , σ k 2 ∑ j = 1 K π j ⋅ ∏ i = 1 n r N y i r | β j T X i r , σ j 2 .

The M-step yields π ̂ k = ∑ r = 1 R τ k r R , β ̂ k = ∑ r = 1 R τ k r ∑ i = 1 n r w i X i r X i r T − 1 ∑ r = 1 R τ k r ∑ i = 1 n r w i y i r X i r , σ ̂ k 2 = ∑ r = 1 R τ k r ∑ i = 1 n r w i ⋅ y i r − β ̂ k T X i , r 2 ∑ r = 1 R n r ⋅ τ k r .

Next, for prediction, on one hand, if the region of a new data point X ̃ i is not known, we estimate the log intensity y ̂ i | X ̃ i = ∑ k = 1 K π k ⋅ X ̃ i T β k . On the other hand, if its region is r _i , then y ̂ i | r i , X ̃ i = ∑ k = 1 K τ r i , k ⋅ X ̃ i T β k .

The mixture model MM-R given in Section 3.3 requires all the flare events of an active region to follow the same regression pattern. This condition can be too restrictive. Note that in our dataset, each flare occurs at a different location and time within an active region, and we can think of the entire data being a collection of events from many active regions. Now, it is reasonable to accommodate the possibility that the heterogeneous nature extends further into each individual flare event within an active region. To model this behavior, we can assign a latent variable z i r to each data point i but impose that z i r ∼ Cat ( ⋅ | π i r ) , where π r ∈ R K and ∑ k = 1 K π k r = 1 . The parameter π k r captures a regional inclination for certain categories of the flaring mechanism. If π ^r (s) are extreme, e.g., taking value 1 in one coordinate and zeros elsewhere, the model MM-H is reduced to the mixture model MM-R given in the previous section. Note that this type of model is sometimes called a mixed membership model in the statistical learning literature (Airoldi et al., 2014), where an active region is a group of members, where the members, in this case, are its flare events. As discussed previously, the weighted complete log-likelihood is optimized to combat the unbalanced data.

Mathematically, the model is parameterized by β 1 , … , β K , σ 1 2 , … , σ K 2 , π 1 , … , π R , where K is the number of global parameters and R is the number of active regions. For active region r = 1, … , R, z i r ∼ Cat ⋅ | π r , i = 1 , … , n r y i r | z i r = k , X i r ∼ N ⋅ | β k T X i r , σ k 2 .

Figure 3 is the probabilistic graphical model representation of model MM-H. There are 2K “global” parameters { β k , σ k 2 } k = 1 K . With the plate notation, there are now n = ∑ r = 1 R n r latent variables z i r for each flare event ( X i r , y i r ) of active region r. Under each active region r, the latent variable z i r is a discrete random variable which takes values in {1, … , K} with probability weight π ^r . Compared to the model described in Section 3.3, MM-R only has R latent variables z ^r , with one “global” weight π. MM-H provides each active region the flexibility of having its own categorical weight π ^r over K linear mechanisms.

Similar to the previously discussed model, the likelihood is p ( y i r | x i r ; π , β , σ 2 ) = ∑ k = 1 K π k r ⋅ N ⋅ | β k T x i r , σ k 2 , and the weighted expected complete log-likelihood optimization problem is equivalent to argmax β , π , σ 2 ∑ k = 1 K ∑ r = 1 R ∑ i = 1 n r E z i r | y i r , X i r ⋅ log ⁡ π k r − w i 2 log ⁡ σ k 2 − w i 2 σ k 2 ⋅ y i r − β k T X i r 2 .

Under these specifications, for parameter estimation, the iterative E-step computes τ i , k r = P z i r = k | y i r , X i r = π k r ⋅ N y i r | β k T X i r , σ k 2 ∑ j = 1 K π j r ⋅ N y i r | β j T X i r , σ j 2 , while the M-step performs the updates π ̂ k r = ∑ i = 1 n r τ i , k r n r β ̂ k = ∑ r = 1 R ∑ i = 1 n r τ i , k r ⋅ w i ⋅ X i r X i r T − 1 ∑ r = 1 R ∑ i = 1 n r τ i , k r ⋅ y i r X i r σ ̂ k 2 = ∑ r = 1 R ∑ i = 1 n r τ i , k r w i ⋅ y i r − β ̂ k T X i , r 2 ∑ r = 1 R ∑ i = 1 n r τ i , k r .

Finally, to perform prediction for a new data point X ̃ i of the region r _i , given r i , X i ̃ , we take y ̂ i ≔ ∑ k = 1 K π k r ⋅ β k T . X i ̃ .

Note that when mixture component K = 1, both models given in Section 3.3 and Section 3.4 reduce to a weighted linear regression model. To perform model selection, the test set is fixed, and the original train set is randomly split into a sub-train set and a validation set for 100 repetitions. For each repetition, we train weighted linear regression (Section 3.1), MM-R (Section 3.3), and MM-H (Section 3.4) on the sub-training set and then apply them to the validation set to obtain the box plots given in Figure 4. In the figures, each box plot visualizes the minimum, first quartile, median, third quartile, and maximum of the collection of generated metrics over 100 repetitions. The averages are also depicted as green dots. The numbers in the x-axis are the number of mixture components to be selected. For MM-R, setting K = 3 yields the best f1 score. A similar conclusion for K can be observed from the breakdown of the rmse for B/C and M/X categories.

The B/C rmse is at the lowest for K = 3, 4, 5, and the M/X rmse is at the lowest for K = 3, 4, as shown in Figures 4C, D. Taking all observations into consideration, we pick K = 3 for MM-R. Following a similar reasoning, we pick K = 3 for MM-H. Note that we combine the rmse of M/X categories and B/C categories since the M/X represents strong flares and B/C represents weak flares. In addition, we are most interested in early warnings of strong flares.

In this section, we demonstrate the clustering results of the mixture model MM-R discussed in Section 3.3 and MM-H discussed in Section 3.4 after training on the observed data. We inspect clustering effects on the marginal space of response y and the marginal space of covariate X, and then we explore what clustering structure entails the interactions between y and X.

Note again that when the number of mixture components K = 1, the models discussed in Section 3.3 and Section 3.4 are just the weighted linear regression model. To produce the figures given in Tables 3–6, we train standard linear regression, weighted linear regression, and models MM-R (Section 3.3) and MM-H (Section 3.4) on the training data, as described in the data section. We run MM-R and MM-H for 100 replications, under each of which training and testing datasets are randomly split, as described in Section 2. We then take the averages as the final results. For each replication, we run our EM procedures five times and select the one with the highest likelihood as the fitted model. We use the validation set to pick the best number of components K. Even though the models predict continuous responses, and we can assess their rmse, to help illustrate the model performances against data imbalance, of which the rmse is not particularly helpful, we compute the f1 score. Since this is a classification metric, continuous responses need to be converted into binary outcomes. We use the same validation set to pick the best binary splitting thresholds. The estimated binary responses are compared with the true labels of strong flares (M/X) and weak flares (B/C) in the data.

The numbers given in Tables 3–6 show that mixture models MM-H (Section 3.4) and MM-R perform similarly, although the former is marginally better. MM-R and MM-H (Section 3.4) significantly outperform the weighted linear regression. This result implies that adding more components improves the performance and, thus, supports the heterogeneous nature of the data. Model MM-H offers more flexibility by extending the heterogeneity pattern to individual flare events. However, interestingly, this flexibility does not noticeably improve the defined metrics. This suggests that the heterogeneity signal is most noticeable at the active region level. We also observe that the performance degrades as the predicting time windows increase, which is expected. Finally, we remark that our predictive performance in the f1 metric (0.383) for 6-h data is lower than that observed by Chen et al. (2019). This is not surprising as we assume linear relationships between covariates and responses and do not account for the temporal nature of solar flares, where future events may correlate with past events. In contrast, Chen et al. (2019) constructed a sophisticated LSTM neural network to extract complex nonlinear signals from the data. We describe potential future directions for improving model performance in the last section. However, the main contribution of this work is demonstrating how mixture models can cluster solar active regions based on the interaction mechanisms between their SHARP covariates and corresponding intensity responses, thus characterizing the heterogeneous nature of active regions.

Both models MM-R 3.3 and MM-H 3.4 perform similarly regarding f1 metrics, even though MM-H is marginally better. As mentioned in the last section, this result seems to suggest that flares from the same AR are intrinsically homogeneous in nature or homologous, as known to the solar physics community (Manchester, 2003; Sui et al., 2004; Liu et al., 2014; Romano et al., 2018). In contrast, individual ARs are most often heterogeneous. As such, we use MM-R for the case studies in this section. As stated in the introduction, the goal of this project is to characterize the heterogeneity among ARs. The cluster membership of the fitted models is particularly interesting, because the model MM-R groups similar active regions together in each cluster.

Because the training and testing sets are created randomly, mixture memberships might change in different replications of the model fitting procedure. Thus, to analyze the cluster assignment in a robust manner, we run the model fitting procedure for MM-R for 100 repetitions. In each repetition, active regions are allocated to different mixture clusters. To standardize the meaning of clusters across repetitions, we assign labels “H,” “L,” and “I” to mixture clusters of which the first quartile of its collections of log intensities is greater than −5.75, smaller than −6.0, and between [−6, −5.75], respectively. By the design of the log intensity threshold for “H,” “I,” and “L,” active regions allocated to “H” labels should be reasonably active in terms of strong flare events, active regions to allocated “L” should be relatively quiet, and “I” in between.

We observe that some active regions have the same labels for all 100 repetitions. For example, ARs 11,124 and 11,109 are assigned to the label “L” in each of the 100 iterations. On the other hand, other active regions might be allocated to different labels over 100 repetitions. For instance, AR 11,967 is assigned to cluster “H” 87% and “I” 13% of the 100 repetitions or AR 12192 85% “H” and 15% “I.” The reason why an active region may have different labels for different repetitions has to do with the randomness of data splitting. In particular, AR 11153 has 28 flare events, of which only one is M class and the rest are B and C ones. As such, if the M-class event is included in the training, it will skew the average log intensities higher than otherwise, and so, it affects its mixture cluster assignment.

The complete list of active region membership is provided in Supplementary Material. Here, we mention the top two ARs for each of the labels “H,” “I,” and “L” in terms of the highest number of label assignments in 100 repetitions. Regarding the “L” label, ARs 11,117 and 11,109 are assigned to “L” 99% and 100% of the times, respectively. For “H,” ARs 11,967 and 12,242 are assigned to “H” 87% and 88% of the times, respectively. For “I,” ARs 11,261 and 11,087 are allocated to “I” 62% and 59% of the 100 repetitions, respectively. Some existing works in the space weather literature corroborate that ARs 11,967 and 12,242 were known to produce strong flares (Solovev et al., 2019; Durán et al., 2020; Joshi et al., 2021). A common trait of ARs with the “I” label is that they have few strong flares among a majority of quiet flare events. In contrast, ARs with “L” labels have only quiet flare events.

To complete our case study, we provide the HMI image and the temporal evolution of each SHARP parameter for the strongest and weakest flare events for AR 11,967 during its existence from January 27 2014 to February 09 2014 in Figure 16. Appendix B provides the same plots for all ARs 11,117, 11,261, 11,967, and 12,242.

a) The left two images are radial components, and the right two ones are horizontal field components in AR 11967 (H label) at its strongest (2014.01.30_16:11:00) and weakest (2014.01.30_13:36:00) flare events during its existence. Arrows show the direction and relative magnitude of the horizontal magnetic field component.

b) (Left) Evolution of SHARP parameters during the strongest (red) and weakest (green) flare events in AR 11967 (H label) and (Right) MM-R estimated β ̂ k values. The dashed lines are the times at their peak intensities. The solid lines are the 6-hour-before for covariates’ regression. The green color is the weakest flare, and red is the strongest one.

All flaring active regions share common basic features, which is a nonpotential magnetic field forming a filament channel over a well-defined polarity inversion line (Green et al., 2018). Beyond this basic feature, there are many possible avenues to eruption. Here, we summarize the observed magnetic structure and evolution of cluster members to determine whether there are features or processes responsible for the heterogeneity of the mixture models. The H cluster contains numerous X-class flares, which received considerable attention in the published literature. The conditions of AR 12,242 leading to the X1.8 flare on December 20 2014 are particularly well described. For example, Solovev et al. (2019) describes the convergence of magnetic flux toward the polarity inversion line leading to both a local and a total maximum gradient of the magnetic field at the time of the flare. Solovev et al. (2019) further described the formation of a magnetic flux rope by reconnection between the converging/colliding sunspots, which erupts to produce the flare. This flaring process has been developed by numerous authors, e.g., Chintzoglou et al. (2019) and Liu (2020). In the case of AR 11,967, a series of flares occurred at a sunspot light bridge, a region of extremely intense and highly sheared magnetic fields produced by flux emergence (Kawabata et al., 2017; Durán et al., 2020). In both examples, we find an intensification of the magnetic field.

For the I cluster, we again find well-described events showing consistent patterns of evolution leading to the flares. The strongest flare from AR 11,087, a C2.7 event, occurred on July 13 2010 at 10:51:00 UT, which is described by Joshi et al. (2015). The flare is characterized by the activation and partial eruption of an active region filament that produces a pair of flair ribbons. Another member of the I cluster, AR 11,261, produced a series of flares, including 4 M-class events, which arose from a complex system of 4 sunspots, one being in a delta configuration (Thalmann et al., 2016). Ye et al. (2018) found that shearing of the photospheric magnetic field associated with flux emergence was a key driver of flares. These observations are consistent with Lorentz force-driven shear flows powering solar eruptions (Ward, 2001; Manchester, 2003; Manchester et al., 2004; Ward, 2007; Fang et al., 2010). Similarly, Sarkar et al. (2019) found that the free energy necessary for flares from active region 11,261 was provided by the shearing motion of moving magnetic features of opposite polarities near the polarity inversion line. The authors also found patterns in the time evolution of the net Lorentz force associated with solar flares.

As stated earlier, the L cluster is dominated by weak flares, which poses two challenges. First, the events are so low in energy that they often occur without significant changes in the photospheric magnetic field and without clear precursors. Second, these low-energy events are far less documented in the literature. However, active region 11,117 is an exception being described in detail in a series of papers by Jiang et al. (2012); Jiang et al. (2016); and Jiang et al. (2017), which we recount here. This region produced a series of small B-class flares observed on October 25 2010 (a date in between our strong and weak flare events), culminating in a C2.3-class event. As described by Jiang et al. (2012), the coronal loops of active region 11,117 (observed by AIA-171) remained largely unchanged, but a flare reconnection was observed at the location of a magnetic null derived from their nonlinear force-free field (NLFFF) model. This eruption event is consistent with topologically driven reconnection models (Titov et al., 2010; Liu et al., 2016). In this case, observations show shear and rotational motions at the photosphere, providing a clear buildup of energy preceding the flares, but no sudden changes precipitate a flare.