Image patch analysis of sunspots and active regions. I. Intrinsic dimension and correlation analysis

Complexity of an active region is related to its flare-productivity. Mount Wilson or McIntosh sunspot classifications measure such complexity but in a categorical way, and may therefore not use all the information present in the observations. Moreover, such categorical schemes hinder a systematic study of an active region's evolution for example. We propose fine-scale quantitative descriptors for an active region's complexity and relate them to the Mount Wilson classification. We analyze the local correlation structure within continuum and magnetogram data, as well as the cross-correlation between continuum and magnetogram data. We compute the intrinsic dimension, partial correlation, and canonical correlation analysis (CCA) of image patches of continuum and magnetogram active region images taken from the SOHO-MDI instrument. We use masks of sunspots derived from continuum as well as larger masks of magnetic active regions derived from the magnetogram to analyze separately the core part of an active region from its surrounding part. We find the relationship between complexity of an active region as measured by Mount Wilson and the intrinsic dimension of its image patches. Partial correlation patterns exhibit approximately a third-order Markov structure. CCA reveals different patterns of correlation between continuum and magnetogram within the sunspots and in the region surrounding the sunspots. These results also pave the way for patch-based dictionary learning with a view towards automatic clustering of active regions.


Introduction
Active regions (AR) in the solar atmosphere have intense and intricate magnetic fields that emerge from subsurface layers to form loops which extend into the corona. When active regions undergo external forcing such as flux emergence and rearrangement, the system may destabilize. The stored magnetic energy is then suddenly released as accelerated particles (electrons, protons, ions) and an increase in radiation called a flare is observed accross the entire electromagnetic spectrum (Phillips, 1991).
The morphology of sunspots is correlated with flare occurence and has therefore received a lot of attention. The Mount Wilson classification scheme (Hale et al., 1919) groups sunspots into four main classes based on the magnetic structure, that is, on the relative locations and sizes of concentrations of opposite polarity magnetic flux. The sunspots with simplest morphology belong to the unipolar α and the bipolar β groups. More complex morphologies are described as βγ when a bipolar sunspot is such that a single north-south polaritiy inversion line cannot divide the two polarities. When a βγ sunspot group contains in addition a δ spot, that is, umbrae of different polarities inside a single penumbra, it is labelled as a βγδ group. The presence of a δ configuration, where large values of opposite polarity exist close together, was identified as a warning of the build up of magnetic energy stress with an increased probability of a large flare (Mayfield and Lawrence, 1985;Sammis et al., 2000). McIntosh (1990) proposes another classification scheme containing 60 classes, thus describing the magnetic structure in greater details. The McIntosh classification is the basis for several flare forecasting methods which estimate the flare occurence rate from historical records of flares and active region classes (Bornmann and Shaw, 1994), possibly combining such information with observed waiting time distribution between flares (Gallagher et al., 2002;Bloomfield et al., 2012).
The McIntosh and Mount Wilson classifications are in general carried out manually, and this results in inconsistencies that stem from human observation bias as well as non-reproducible catalogs. To overcome these caveats, some supervised machine learning methods have been proposed to automatically classify sunspot groups according to these schemes. Stenning et al. (2013) extract various measurements from continuum and magnetogram images, and then feed these into a machine learning classifier which reproduces the Mount Wilson classification. Colak and Qahwaji (2008) employ neural networks and supervised classification techniques to reproduce the McIntosh scheme and use those results in a flare forecasting system (Colak and Qahwaji, 2009). While these approaches reduce the human bias, they do not use the information present in sunspot images in an optimal way and make the study of AR dynamic behavior impractical.
Several attempts were made to find quantitative descriptors of an active region's complexity. McAteer et al. (2005) showed that fractal dimension of an active region alone cannot distinguish between the various Mount Wilson classes. The generalization to multifractal spectrum, where each scale has its own fractal dimension, allowed to study in greater details the evolution of active region in view of distinguishing between quiet and flare-productive active regions. Box counting method (Georgoulis, 2005;Abramenko, 2005;Conlon et al., 2008) as well as more acurate methods based on continuous wavelet transform (Kestener et al., 2010;Conlon et al., 2010) were employed. Continuous wavelet transforms and energy spectrum were also used with a similar purpose in Hewett et al. (2008); . Wavelet basis functions act as a microscope to describe local discontinuities and gradients in an image, and Ireland et al. (2008) used two multiresolution analyses to compute at various length scales the gradients of the magnetic field along lines separating opposite polarities. Using a dataset of about 10 000 magnetogram images, they showed that, at all length scales, those gradients increase going from α to β, βγ, and βγδ classes.
However, a wavelet analysis is known to generate artifacts due to the particular shape of the specific wavelet functions. To remedy this problem, dictionary learning methods based on patch analysis have been proven to be more flexible and provide superior results for problems such as, e.g., image denoising (Elad and Aharon, 2006). Technically, a patch is a m × m-pixel neighborhood, and a patch analysis of a n-pixel image will process the m 2 × n data matrix that collects the overlapping patches. Figure 1 presents such a patch and its column representation.
In this paper, we carry out a patch analysis of a set of sunspots and active region magnetogram images that span the four main Mount Wilson classes. We estimate the intrinsic dimension of the local patches, and show how it relates to the Mount Wilson classification. We also study patterns of local correlation using partial correlation and canonical correlation analysis, which reveal some characteristics of simple and more complex active regions. Such analysis also serves as a preparation to an unsupervised clustering of active region using patch-based dictionary learning which will be presented in a companion paper.
Section 2 describes our dataset. Unlike previous works, our approach combines information from two modalities: photospheric continuum images and magnetograms, both obtained by the Michelson Doppler Imager (MDI) onboard the Solar and Heliospheric Observatory (SOHO) spacecraft. We consider 424 active regions spanning the four main Mount Wilson classes. We use SMART masks (Higgins et al., 2011) to delineate the boundaries of magnetic active regions, and the STARA algorithm (Watson et al., 2011) which provides masks for umbrae and penumbrae from the continuum images. These two masks enable us to differentiate between pixels belonging to the actual sunspots and pixels featuring the region surronding the sunspots.
In Section 3, the intrinsic dimension of the image patches extracted from the two modalities is estimated using both linear and non-linear methods. The linear method relies on Principal Component Analysis (PCA) (Jolliffe, 2002), while the non-linear method relies on a k-Nearest Neighbor graph approach (Costa and Hero III, 2006;Carter et al., 2010). We show that the intrinsic dimension is related to the complexity of the sunspot groups. Section 4 identifies the spatial and modal interactions of the patches at different scales by estimating the partial correlation and by using canonical correlation analysis (CCA) (Muller, 1982;Nimon et al., 2010). This gives insight about relationships that may exist between active region complexity and the correlation patterns.
This paper expands and refines some of the work in . Whereas  used fixed size square pixel regions centered on the sunspot group as input to the analyses, in this paper SMART detection masks are used. A larger set of images is considered in all methods which enables us to analyze the relationships of intrinsic dimension and correlation with AR complexity. We also explore the partial correlation of patches which was not included in .

Data
The data used in this study are taken from the Michelson Doppler Imager (MDI) instrument (Scherrer et al., 1995) on board the SOHO Spacecraft.
Within the time range of 1996-2010, we select a set of 424 ARs as follows. Using the information from the Solar Region Summary reports compiled by the Space Weather Prediction Center of NOAA http://www.swpc.noaa.gov/ftpdir/forecasts/SRS/, we consider ARs located within 30 • of the solar meridian. We looked at a maximum of two-hundred instances per Mount Wilson types α, β, βγ, and βγδ. Out of this first selection, we removed AR with a longitudinal extent smaller than four, and finally we checked if MDI continuum and magnetogram data were available. This provides us with a number of ARs in each Mount Wilson class as given by Table 1. In our analysis, we also divide the ARs into two groups: simple ARs (α and β) and complex ARs (βγ and βγδ).
AR are observed using two modalities: photospheric continuum images and magnetogram. SOHO-MDI provides almost continuous observations of the Sun in the white-light continuum, in the vicinity of the Ni i 676.78 nm photospheric absorption line. These photospheric intensity images are primarily used for sunspot observations. MDI data are available in several processed "levels". We use level-1.8 images, and rotate them with North up. SOHO provides two to four MDI photospheric intensity images per day with continuous coverage since 1995. Using the same instrument level, 1.8 line-of-sight (LOS) MDI magnetograms are recorded with a nominal cadence of 96 minutes. The magnetograms show the magnetic fields of the solar photosphere, with negative (represented as black) and positive (as white) areas indicating opposite LOS magnetic-field orientations.
As stated in Section 1, SMART masks (Higgins et al., 2011) are used to determine the boundaries of magnetic active regions from MDI magnetograms. Those masks are applied also on continuum images to determine the surrounding part of the sunspot that is affected by magnetic fragments as seen in magnetogram images. Similarly, the STARA algorithm (Watson et al., 2011) provides masks for sunspots (umbrae and penumbrae) from MDI continuum and those masks are applied on magnetogram images to determine the AR cores corresponding to the sunspots. Combining these two types of masks provides thus two sets of pixels within each AR: those belonging to the sunspots themselves as found by STARA and those belonging to the magnetic fragments (or background) within an AR as found by the difference set between the SMART and STARA masks.
As in  we use image patch features to account for spatial dependencies using square patches of pixels. Thus if a SMART mask of an image has n pixels and we use a m × m patch, the corresponding continuum data matrix X is m 2 × n where the ith column contains the pixels in the patch centered at the ith pixel. The magnetogram data matrix Y is formed in the same way and the full data matrix is Z = X Y with size 2m 2 × n. We let z i denote the ith column of Z. The images from both modalities are also normalized prior to analyzing them. Since some sunspot and active region features can be quite small and to limit the effects of high dimensionality on our analysis, we primarily use 3 × 3 patches in each modality although larger patches are used in Section 4 when analyzing spatial correlations in the images.

Intrinsic Dimension Estimation
The goal of this section is to determine the number of intrinsic parameters or degrees of freedom required to describe the spatial and modal dependencies using image patches. We consider 3 × 3 patches within both the continuum and magnetogram images giving an extrinsic dimension of 18. It is also important to know whether linear analyses can be accurately applied to the data or whether non-linear techniques are required. To answer this question, we estimate the local intrinsic dimension using a method appropriate for linear subspaces and a method appropriate for any (linear or non-linear) smooth subspace and then compare the results.

PCA: A Linear Estimator
Principal Component Analysis (PCA) (Jolliffe, 2002) finds a set of linearly uncorrelated vectors (principal components) that can be used to represent the data. PCA has been used previously for various purposes in solar-physics and space-weather literature, e.g. to study the background and sunspot magnetic fields (Lawrence et al., 2004;Cadavid et al., 2008;Zharkova et al., 2012), for analysis of solar wind data (Holappa et al., 2014), or to reduce dimensionality (Dudok DeWit and Auchère, 2007).
In PCA, the principal components are the eigenvectors of the covariance matrix Σ: where x and y are random vectors of dimension 9, x being a patch from the continuum image, and y the corresponding patch from the magnetogram image. The eigenvalues indicate the amount of variance accounted for by the corresponding principal component. A linear estimate of intrinsic dimension is the number of principal components that are required to explain a certain percentage of the variance. By nature, PCA is a global operation and so it provides a global estimate of the intrinsic dimension. We can obtain more local estimates by performing PCA separately on the areas within the sunspots and on the magnetic fragments. These areas are separated using the STARA and SMART masks.

k-NN: A General Estimator
The general method we use is a k-Nearest-Neighbor (k-NN) graph approach with neighborhood smoothing (Costa and Hero III, 2006;Carter et al., 2010). The intuition behind the method is that we grow the k-NN graph from a point z i by adding an edge from z i to z j if z j is within the k nearest neighbors of z i . The growth rate of the total edge length of the graph is related to the intrinsic dimension of the data in a way that enables us to estimate it.
One advantage of the k-NN method is it enables us to estimate the local intrinsic dimension by limiting the growth of the graph to a smaller neighborhood. This provides an estimate of intrinsic dimension at each pixel location in the image which allows us to more easily visualize the intrinsic dimension estimates. Technical explanation of the k-NN method and more details on local intrinsic dimension are given in Appendix A.1.

General Results
Given 3 × 3 patches with extrinsic dimension 18 when using both the continuum and magnetogram images, we estimate the intrinsic dimension within the sunspots and magnetic fragments for all 424 ARs using both the k-NN approach and PCA. Figure 2 shows two examples of the estimated local intrinsic dimension using the k-NN method and the corresponding continuum and magnetogram images. One set of images corresponds to an α group while the other set is a βγδ group. In these examples, areas with more spatial structure, such as within the sunspots, have lower intrinsic dimension. Fewer parameters are required to accurately represent structured data than noise and so the intrinsic dimension is lower. Table 2 provides the mean and standard deviation of the intrinsic dimension estimates within the sunspots and magnetic fragments. These statistics are also provided for ARs within the main Mount Wilson classes (α, β, βγ, and βγδ). We provide PCA results for the cases where we estimate the intrinsic dimension as the number of components required to explain 97% and 98% of the variance, respectively. For the k-NN method, we provide the results in two ways. For one, we take the mean of local intrinsic dimensions within each image (separating the 'sunspot' from the 'magnetic fragments') and then calculate the mean and standard deviation of these means. The statistics in this category correspond to the mean and standard deviation of the average intrinsic dimension of each image and are more directly comparable to the PCA results. However these results may be affected slightly by small sunspot groups. For the other approach, we pool all of the local estimates (again separating sunspots from magnetic fragments) and then calculate the mean and standard deviation. These results correspond to the mean and standard deviation of the pixels within each region and category and are less affected by small sunspot groups. α β βγ βγδ All Sunspots k-NN, pooled 3.9 ± 1.2 4.4 ± 1.0 4.4 ± 0.9 4.5 ± 0.7 4.3 ± 1.0 Sunspots k-NN, means 4.0 ± 1.0 4.8 ± 1.2 4.4 ± 0.6 4.4 ± 0.4 4.5 ± 1.0 Sunspots PCA 97% 3.7 ± 0.8 4.4 ± 0.9 4.3 ± 0.6 4.2 ± 0.5 4.3 ± 0.8 Sunspots PCA 98% 4.5 ± 0.9 5.2 ± 1.0 5.1 ± 0.8 5.0 ± 0.5 5.0 ± 0.9 Fragments k-NN, pooled 8.0 ± 1.1 8.0 ± 1.2 7.7 ± 1.2 7.6 ± 1.2 8.0 ± 1.2 Fragments k-NN, means 8.0 ± 0.4 7.9 ± 0.5 7.6 ± 0.5 7.6 ± 0.5 7.8 ± 0.5 Fragments PCA 97% 7.7 ± 1.6 7.1 ± 1.2 6.2 ± 1.2 6.2 ± 1.3 6.8 ± 1.4 Fragments PCA 98% 9.1 ± 1.6 8.5 ± 1.3 7.5 ± 1.2 7.4 ± 1.3 8.1 ± 1.4 Table 2. Estimated intrinsic dimension results for different groups of ARs in the form of mean±standard deviation. The complex ARs have higher intrinsic dimension within the sunspots than the simple ARs but lower intrinsic dimension within the magnetic fragments.
From Table 2, it is clear that the intrinsic dimension is lower within the sunspots than in the magnetic fragments for all methods. This is expected as there is more spatial structure within the images inside the sunspots than in the magnetic fragments, especially in the continuum image.
The average PCA estimate with a 97% threshold and the average mean k-NN estimate give similar results inside the sunspot while the average 98% PCA estimate is closest to the average mean k-NN estimate within the magnetic fragments. If linear methods were not sufficient to represent the spatial and modal dependencies, we would expect the PCA results to be much higher than the k-NN results when using comparable thresholds as more linear than nonlinear components would be required to accurately represent the data. However, this close agreement between the general and linear results suggests that linear methods are sufficient and that linear dictionary methods would be appropriate for these data.

Patterns Within the Mount Wilson Groups
For both the PCA and k-NN methods, the average estimated intrinsic dimension is lower within the sunspots in α groups than in the more complex groups such as βγδ. This is consistent with Figure 2 and may be related to the lower complexity of α groups. These exhibit more spatially coherent images, which can be described using a lower number of basis elements, and hence have a lower intrinsic dimension.
Within the magnetic fragments, the opposite trend occurs where the less complex groups have higher intrinsic dimension. This suggests that the magnetic fragments are fewer, weaker, and less structured outside of the α and β groups compared to the more complex regions, leading to a more noiselike background in their magnetic fragments. This hypothesis is supported by the normalized histograms of the mean k-NN estimates of intrinsic dimension and the normalized histograms of the pooled k-NN estimates in Figures 3 and 4. The histograms of mean intrinsic dimension show that within the magnetic fragments, α groups generally have higher mean intrinsic dimension than βγδ groups. In fact, no α groups have a mean intrinsic dimension less than 7.5 within the magnetic fragments. However, the normalized histograms of the individual patch estimates show a significant number of patches with intrinsic dimension less than 7.5 within the fragments. This suggests that for each α group, the majority of the patches have higher intrinsic dimension in the magnetic fragments and are thus more noiselike. In contrast, there are some βγδ groups where the mean intrinsic dimension of the mangetic fragments is lower (less than 7.5) and so these magnetic fragments are dominated by patches with more structure.
Table 2 also shows that the standard deviation of the estimates within the sunspots decreases as the complexity increases as measured by the Mount Wilson classification scheme. The histograms in Figures 3 and 4 can be used to determine the cause. From the histograms, it is clear that within the sunspots the intrinsic dimension of α groups does not have a Gaussian distribution. In this case, most of the estimates are between 3 and 5. However, there are a significant number of outliers with intrinsic dimension greater than 5. The presence of these outliers contributes to the high standard deviation. This is in contrast to the intrinsic dimension of βγ and βγδ groups inside the sunspot which have fewer outliers and thus smaller standard deviations.
The outliers in the α groups correspond to small sunspots. The number of pixels within the α sunspots with average intrinsic dimension ≥ 6 range between 10 and 53 with a median of 16. In these cases, the spatial structure of the sunspots may be more similar to the magnetic fragments than the spatial structure of larger sunspots. Thus the intrinsic dimension is higher in the small sunspots.
A similar phenomenon occurs within the β groups. Note that in Table 2, the average and standard deviation of the mean intrinsic dimension of the β groups within the sunspots is higher than for all other groups. This is also caused by a few outliers that have high average intrinsic dimension due to the small size of the sunspots. When individual local intrinsic dimension estimates of the patches Figure 3. Normalized histograms of mean estimated intrinsic dimension of α, β, βγ, and βγδ groups using the k-NN method. The distributions of intrinsic dimension differ by complexity with simpler AR groups having higher (resp. lower) intrinsic dimension within the sunspot (resp. magnetic fragments).
from these small sunspots are pooled with the estimates from all other β patches, the average intrinsic dimension is more aligned with that of the other Mount Wilson types. Additionally, ignoring the biggest outliers in the mean intrinsic dimension (defined as having mean intrinsic dimension > 6.25) gives an average mean intrinsic dimension of 4.6 for the β groups which is more aligned with the other groups.
The distribution of intrinsic dimension within the magnetic fragments also differs by complexity based on Figures 3 and 4. The complex ARs have more patches and images with lower intrinsic dimension than the simple sunspots which is consistent with Table 2.
In summary, based on the estimated intrinsic dimension of the image patches, relatively few parameters are required to accurately represent the data. We have found that the distribution of local intrinsic dimension varies based on the complexity of the sunspot group with the more complex sunspots having higher (resp. lower) intrinsic dimension within the sunspot (resp. magnetic fragments). Additionally, the standard deviation of the intrinsic dimension is higher within the sunspot in the simpler sunspots than the complex ARs. This is due to the presence of small sunspots among the simpler ARs that tend to have less spatial structure and thus a higher intrinsic dimension than typical sunspots. We have also shown that linear methods should be sufficient to accurately analyze the data. . Normalized histograms of pooled local estimates of intrinsic dimension of α, β, βγ, and βγδ groups using the k-NN method. The distributions of intrinsic dimension differ by complexity with simpler AR groups having higher (resp. lower) intrinsic dimension within the sunspot (resp. magnetic fragments).

Spatial and Modal Correlations
As the results in the previous section indicate that linear methods are likely sufficient to represent the spatial and modal dependencies within a sunspot, we analyze the linear correlation over patches. There are several purposes for this analysis. One is we are interested in the relationship between the correlations and sunspot group complexity. A second purpose is to receive some insight into how large of a patch we should use to sufficiently capture the spatial and modal correlations for future analysis. Finally, we would like to analyze the modal dependencies. If the two modalities are independent, then there is no benefit to processing them together. On the other hand, if the two modalities are completely dependent, then we only need to process one of the modalities since the other modality would not contain additional information. We accomplish these purposes by estimating the partial correlation and by performing canonical correlation analysis (CCA) within the sunspots and within the magnetic fragments.

Partial Correlation: Methodology
The partial correlation measures the correlation between two random variables while conditioning on the remaining random variables. The intuition behind partial correlation can be best explained with the linear regression concept. Suppose you want to compute the partial correlation between two Figure 5. Estimated partial correlation matrices of patch data from within the sunspots and the magnetic fragments using 3 × 3 (left) and 5 × 5 (right) patches. The theoretical thresholds (Hero and Rajaratnam, 2011) for significance to attain a 0.05 false alarm rate are 0.0070 and 0.0014 for within the sunspots and magnetic fragments, respectively when using a 3 × 3 patch. For the 5 × 5 patch, the thresholds are 0.0080 and 0.0016, respectively. Statistically insignificant values are set to zero.
variables X 1 and X 2 given a set of variables X. First, compute the linear regression using variables in X to explain X 1 and obtain the associated residuals r X 1 . Proceed similarly for X 2 and get residuals r X 2 . The partial correlation between X 1 and X 2 is then equal to the (usual) correlation between r X 1 and r X 2 , for which the effect of variables X have been removed.
In our context, let x be a patch from the continuum image, and |y| be the magnitude (entrywise absolute value) of the corresponding patch from the magnetogram. The partial correlation matrix P = P xx P x|y| P |y|x P |y||y| and its off-diagonal elements can be derived from the inverse correlation matrix (see Appendix A.2). We use the magnitude of the magnetogram data since both positive and negative polarities affect the continuum image in similar ways. Figure 5 gives the estimated partial correlation matrices when using 3 × 3 and 5 × 5 patches. The patches are extracted from all of the active regions and divided using the STARA and SMART masks into sunspots and magnetic fragments as before. The partial correlation of 3 × 3 patches is quite strong within both modalities. Based on a false alarm rate of 0.05, the theoretical thresholds for significance for the partial correlation (Hero and Rajaratnam, 2011) of the 3 × 3 patches are approximately 0.0070 and 0.0014 for within the sunspots and magnetic fragments, respectively. For the 5 × 5 patches, the thresholds are 0.0080 and 0.0016, respectively. Given these thresholds, the partial correlation is statisitically significant for nearly all values within the modalities (P xx and P |y||y| ) using the 3 × 3 patches. The cross-partial correlation when using the signed magnetic field (i.e. P xy and P yx ) is very near zero in both regions (not shown). However, when we take the absolute value of the magnetogram patches, then the magnitude of the cross-partial correlation (P x|y| and P |y|x ) is much higher in both regions suggesting that the correlation between the modalities is significant. The partial correlation is also stronger in magnitude in all cases within the sunspots than within the magnetic fragments.

Partial Correlation: Results
The partial correlation matrices are very structured. In both sunspots and magnetic fragments, the pentadiagonal-like structure within the modalities suggests that the image is generally stationary Figure 6. Partial correlation patches extracted from the columns in the sunspot partial correlation matrix corresponding to the center pixels. The partial correlation is stronger within the continuum. Figure 7. Partial correlation matrices within the sunspots using the data from α (left) and β ARs (center). Statistically insignificant values are again set to zero. (Right) difference between the absolute value of the α and β matrices. The α sunspots are more (resp. less) strongly correlated within (resp. between) the modalities than the β groups.
with approximately a third order nearest neighbor Markov structure. Such structure is clearly seen in the matrices for 5×5 patches. The cross-partial correlation also has a pentadiagonal-like structure although the correlation is not as strong as within the modalities.
To better see the spatial correlations, in Figure 6 we plot the partial correlation patches taken from the columns of the sunspot partial correlation matrix corresponding to the center pixels when using 5 × 5 patches. From these images, the greater partial correlation within the continuum is even more evident. It's also evident that the correlation is slightly higher in magnitude in the vertical direction than the horizontal direction. Nearly all sunspots in this study are located within (−30 • , +30 • ) from both the central meridian and the equator, and so projection effect are unlikely to cause this difference. The difference in correlation may be a feature of the sunspots themselves, but this may be difficult to determine since the difference in partial correlation is small. Some slight differences exist in the partial correlation matrices restricted to certain Mount Wilson classes. As an example, Figure 7 contains the partial correlation matrices within the sunspots after restricting the data to α and β groups as well as the difference between the absolute value of the two matrices. The α partial correlation matrix is higher in magnitude within the modalities than the β matrix but lower between the modalities. Within modalities, the strongest differences (a maximum of 0.056 and 0.067 within the continuum and magnetogram, respectively) are in the entries that correspond to pixels that are farther away from each other. In contrast, within the cross-partial correlation, the strongest differences (a maximum of 0.072) between the two AR types are in the entries that correspond to pixels that are close to each other. A similar pattern holds when comparing the α matrix to the matrices of the more complex groups.
Overall, the partial correlation matrices indicate that no larger than a 5 × 5 patch is necessary to capture the local spatial dependencies. However, in the remainder of our analysis, we choose a 3 × 3 patch for several reasons. One reason is that there are some small sunspots that we include in our data. A patch that is too big would be unable to capture the features of these small sunspots. Additionally, we would like to control the size of the extrinsic dimension of our data (the sum of the total number of pixels in both patches) since this affects the accuracy of the analysis. Given these concerns, we see that 3 × 3 patches capture most of the spatial correlation. This is evident from Figure 5 where the partial correlation between pixels on opposite corners of a 3 × 3 patch is near zero and other pixels that are similarly far away from each other have low partial correlation. Thus a 3 × 3 patch strikes a good balance between scale, extrinsic dimension, and capturing the spatial correlation.

Canonical Correlation Analysis: Methodology
To further investigate the correlation between the modalities, we use canonical correlation analysis (CCA) on the continuum patch x and the magnitude of the magnetogram patch y. CCA finds patterns and correlations between two multivariate data sets (Muller, 1982;Nimon et al., 2010) and was used previously in the context of space weather e.g. for the combined analysis of solar wind and geomagnetic index data sets (Borovsky, 2014).
In our application, CCA finds vectors a i and b i such that the correlation ρ i = corr(a T i x, b T i |y|) is maximized and the pair of random variables u i = a T i x and v i = b T i |y| are uncorrelated with all other pairs u j and v j , j i. The variables u i and v i are called the ith pair of canonical variables while the vectors a i and b i are the canonical vectors. The solution a i is the ith eigenvector of the matrix Σ −1 xx Σ x|y| Σ −1 |y||y| Σ |y|x . The vector b i is found similarly (Härdle and Simar, 2007).

Canonical Correlation Analysis: Results
Here we focus on 3 × 3 patches and apply CCA to all 424 image pairs using the magnitude of the magnetogram patches. Figure 8 (left and center) shows histograms of the estimated values of ρ 1 . Within the sunspots, there are many groups with a near perfect correlation between the modalities and none of the groups have an estimated value below 0.41. The right plot in Figure 8 plots the estimated values of ρ 1 vs. the number of samples used within the sunspots. Based on this plot, there are many ARs with high correlation and few patch samples suggesting that the correlation may be spurious. However, all of the estimated values are statistically significant as defined by the threshold given by Hero and Rajaratnam (2011) using a false alarm rate of 0.05 (shown as the magenta line in Figure 8). The histogram of ρ 1 within the magnetic fragments (Fig. 8, center) is quite different from the sunspot histogram (Fig. 8, left). ρ 1 is generally lower within the magnetic fragments than within the sunspots which is consistent with the results in Figure 5. All of the ρ 1 values are statistically significant.
The distributions of ρ 1 differ slightly when comparing simple sunspot groups (α and β) with complex groups (βγ and βγδ). Figure 9 shows that complex groups generally have lower correlation between the modalities within the sunspots than the simpler groups. The estimated Hellinger distance (see Appendix A.3) between the distributions using the divergence estimator in Moon and Hero III (2014a,b) is 0.22. Based on the central limit theorem of the estimator (Moon and Hero III, Figure 8. Histograms of estimated ρ 1 using CCA for within the sunspots (left) and the magnetic fragments (center) using 3 × 3 patches. Right: Scatter plot of ρ 1 values and the number of samples available for within the sunspots. All points are above the magenta line which gives the threshold for statistical significance at a false alarm rate of 0.05 (Hero and Rajaratnam, 2011). Figure 9. ρ 1 histograms of complex (βγ and βγδ) and simple (α and β) regions within the sunspots (left) and the magnetic fragments (right). The simple ARs are generally more correlated within the sunspots but less correlated within the magnetic fragments. The difference between the sunspot distributions, as measured by the Hellinger distance, is statistically significant. 2014b), this value is statistically significant with a p-value of 1.6×10 −12 . At least some of this difference is likely due to the smaller size of the simpler groups (and thus smaller sample size). However, it is unlikely to fully explain the difference given that there are many simple sunspot groups with high correlation and sufficient sample size.
Within the magnetic fragments, there are many more simple regions than complex regions with ρ 1 < 0.4 (see the histogram in Figure 9, right). This could be related to the same phenomena that causes the intrinsic dimension to be higher within the magnetic fragments of simple sunspot groups observed in Section 3.3. However, the estimated Hellinger distance between these distributions is 0.016. Using the same statistical test, this estimate is not statistically significant with a p-value of 0.31. Thus the distributions are not statistically different from each other.
To analyze the spatial patterns that produce the highest correlation between modalitites, we apply CCA to the entire dataset. Figure 10 plots ρ i for i = 1, . . . , 9 for within the sunspots and within the magnetic fragments. The ρ i are all statistically significant. Notice that the ρ i are higher within the sunspots than the magnetic fragments which is consistent with the results in Figures 5 and 8.  Figure 10 shows the canonical patches a i and b i for i = 1, . . . , 6 when using all the data from within the sunspots. These are the spatial patterns within the two modalities that are most correlated with each other. The canonical patches have a "saddle-like" appearance where the gradient is positive in some directions and negative in others. For example, in a 4 , the pixels to the left and right of the center are very negative but the pixels in the corners are all very positive. Note that these vectors corespond to centered values with respect to the mean patches.
Comparing the a i s to the b i s shows that the b i s are approximately equal to the negative of the a i s. This makes sense as sunspots within the continuum images correspond to a decrease in value relative to the background while ARs within the magnitude of the magnetogram images correspond to an increase in value relative to the background.
We also performed CCA separately on the data from the Mount Wilson classes. Figure 11 plots the ρ i values for each class and the first canonical patches a 1 and b 1 . For ρ 1 and ρ 2 , the values for each class decrease in order of complexity (α, β, βγ, βγδ). This is consistent with our comparison of the partial correlation matrices in Figure 7 where the partial correlation was generally higher (in magnitude) for the α groups than the others. This is also consistent with the intrinsic dimension analysis in Section 3 where the intrinsic dimension generally increases with complexity. This is because if the correlation between and within modalities is higher, then fewer parameters are required to accurately describe the data which results in a lower intrinsic dimension.
The canonical patches a 1 and b 1 have similar patterns across the different classes although the patches for the βγ class are flipped compared to the others. The magnitude of the values in the βγδ patches are also smaller than the those of the other patches.
Overall, the results of this section suggest that the two modalities are correlated in both the sunspots and the magnetic fragments and are therefore not independent. The correlation is stronger within the sunspots compared to the magnetic fragments and stronger within the sunspots in simple ARs compared to complex ARs. However, the correlation is not perfect and so there may be an advantage to including both modalities in the classification of sunspots and flare prediction.  Figure 10 but the patches differ slightly from class to class.

Conclusion
Existing AR categorical classification systems such as the Mount Wilson and McIntosh schemes describe geometrical arrangements of the magnetic field at the largest length scale. In this work, we have focused on the properties of the ARs at fine length scale. We showed that when we analyze the global statistics or attributes of these local properties, we find differences between the simple and complex ARs as defined using the large scale characteristics. So by this approach, we are analyzing both the large and fine scale properties of the images. Such results might be due to the multi-scale properties of the magnetic fields, as evidenced previously in Ireland et al. (2008).
The local intrinsic dimension based on the k-NN approach combines both continuum and magnetogram observations and provides some measure of local regularity for those images. Further differences between the Mount Wilson classes may be found by comparing the histograms or distributions of local intrinsic dimension of each individual AR instead of only comparing the means or pooled estimates as we did in this paper. There are several options to perform such comparisons. Each histogram could be treated as a vector, or we could consider the underlying probability density function within the framework of functional analysis. Supervised (using Mount Wilson classes) or unsupervised classification could be performed. Another option would be to view the set of histograms belonging to a specific class as samples from a distribution of vectors (or a distribution of probability density functions). Different classes could then be compared using divergence measures such as the Hellinger distance described in Appendix A.3. This work also highlighted specific behaviors of the core of active regions (that corresponds to the sunspot masks in continuum) and magnetic fragments (the surrounding part of AR), as well as the difference of these two regions as a function of the Mount Wilson classification. We found that within the sunspots, the spatial and modal correlations are stronger than within the magnetic fragments. Additionally, simpler ARs were found to have higher correlation between the modalities within the sunspots than the complex ARs.
This study paves the way for further analysis based on dictionary learning. Knowledge of the intrinsic dimension allows us to choose the dictionary size. Moreover the results of Section 3 showed that linear dictionary learning methods are sufficient. The spatial and modal correlation analysis in Section 4 justifies a choice of a patch size of 3 × 3 and confirms that both modalities (continuum and magnetogram) should be used in dictionary learning.

A.1. Intrinsic Dimension Estimation of Manifolds
Consider data that are described in an extrinsic Euclidean space of d dimensions. However, suppose the data actually lie on a lower dimensional manifold M. Thus the intrinsic dimension m of the data corresponds to the dimension of M. For example, data may be given to us in a 3 dimensional space but lie on the surface of a sphere. Thus the intrinsic dimension of the data would be 2.
In some cases, data points from the same dataset may lie on different manifolds. For example, part of the data with an extrinsic dimension of 3 could lie on the surface of a sphere (m = 2) while another part may lie on a circle (m = 1). We then say that data points from these different manifolds have a different local intrinsic dimension. The local intrinsic dimension gives some measure of the local complexity of the image. Additionally, the local intrinsic dimension is useful for dictionary learning because we can use it to determine whether different-sized dictionaries should be used for different regions, e.g. within the sunspots and outside of the sunspots.
We now describe the k-NN estimator of intrinsic dimension in more detail. For a set of independently identically distributed random vectors Z n = {z 1 , . . . , z n } with values in a compact subset of R d , the k-nearest neighbors of z i in Z n are the k points in Z n \{z i } closest to z i as measured by the Euclidean distance || · ||. The k-NN graph is then formed by assigning edges between a point in Z n and its k-nearest neighbors. The intrinsic dimension is related to the total edge length of the k-NN graph and can be estimated based on this relationship. The k-NN graph is then formed by assigning edges between a point in Z n and its k-nearest neighbors and has total edge length defined as where γ > 0 is a power weighting constant and N k,i is the set of k nearest neighbors of z i . It has been shown that for large n, L γ,k (Z n ) = n α(m) c + n , where α = (m − γ)/m, c is a constant with respect to α(m), and n is an error term that decreases to zero a.s. as n → ∞ (Costa and Hero III, 2006). A global intrinsic dimension estimatem is found based on this relationship using non-linear least squares over different values of n (Carter et al., 2010). A local estimate of intrinsic dimension at a point z i can be found by running the algorithm over a smaller neighborhood about z i . The variance of this local estimate is then reduced by smoothing via majority voting in a neighborhood of z i (Carter et al., 2010).

A.2. Partial Correlation
Let z be a random vector with size m. Let Σ be the covariance matrix of z, that is Σ i j = Cov(z i , z j ), and let K = Σ −1 be the inverse of the covariance matrix, also called the precision matrix.
The partial correlation between z i and z j given all the other variables z\{z i , z j } measure the degree of correlation between these two variables after removing the effect of the remaining ones. Let P i j denote the partial correlation between z i and z j . It has been shown (Lauritzen, 1996) that P i j can be related to the elements of the precision matrix K as follows:

A.3. Estimating the Hellinger Distance
Information divergences are a class of functionals that measure the difference between two probability distributions. The most popular divergence measure is the Kullback-Leibler divergence (Kullback and Leibler, 1951). The Hellinger distance is another divergence measure and is defined as where f and g are the two probability densities being compared. The Hellinger distance is a metric which is not true of divergences in general. we use the nonparametric divergence estimator derived in Moon and Hero III (2014a,b). In Moon and Hero III (2014a), it was shown that this estimator converges to the true divergence with mean squared error convergence rate O(1/T ) where T is the number of samples from each probability distribution. In Moon and Hero III (2014b), it was shown that the distribution of the normalized version of this estimator converges to the standard normal distribution. We can use this fact combined with a bootstrap estimate (Efron and Tibshirani, 1994) of the variance of the estimator to test the hypothesis that the divergence is zero (and hence the distributions are equal).