Complexity / 2017 / Article

Research Article | Open Access

Volume 2017 |Article ID 3076810 |

Carlos A. L. Pires, Abdel Hannachi, "Independent Subspace Analysis of the Sea Surface Temperature Variability: Non-Gaussian Sources and Sensitivity to Sampling and Dimensionality", Complexity, vol. 2017, Article ID 3076810, 23 pages, 2017.

Independent Subspace Analysis of the Sea Surface Temperature Variability: Non-Gaussian Sources and Sensitivity to Sampling and Dimensionality

Academic Editor: Davide Faranda
Received23 May 2017
Accepted10 Jul 2017
Published22 Aug 2017


We propose an expansion of multivariate time-series data into maximally independent source subspaces. The search is made among rotations of prewhitened data which maximize non-Gaussianity of candidate sources. We use a tensorial invariant approximation of the multivariate negentropy in terms of a linear combination of squared coskewness and cokurtosis. By solving a high-order singular value decomposition problem, we extract the axes associated with most non-Gaussianity. Moreover, an estimate of the Gaussian subspace is provided by the trailing singular vectors. The independent subspaces are obtained through the search of “quasi-independent” components within the estimated non-Gaussian subspace, followed by the identification of groups with significant joint negentropies. Sources result essentially from the coherency of extremes of the data components. The method is then applied to the global sea surface temperature anomalies, equatorward of 65°, after being tested with non-Gaussian surrogates consistent with the data anomalies. The main emerging independent components and subspaces, supposedly generated by independent forcing, include different variability modes, namely, The East-Pacific, the Central Pacific, and the Atlantic Niños, the Atlantic Multidecadal Oscillation, along with the subtropical dipoles in the Indian, South Pacific, and South-Atlantic oceans. Benefits and usefulness of independent subspaces are then discussed.

1. Introduction

The climate system is a highly complex dynamical system with many nonlinearly interacting modes of variability. A characteristic feature of the system is its non-Gaussianity. The nonnormality of climate is an important determinant of extreme events and may provide insights into understanding the system dynamics [13]. For example, it is generally accepted that understanding non-Gaussian statistics of weather and climate has important consequences in atmospheric research not least because weather/climate risk assessment depends on understanding the tail of the system probability density function (PDF). Because of the complex interactions involved in the system, climate signals tend to be mixed. An important problem in climate research is to be able to disentangle main (or independent) signals from the mixture.

Empirical orthogonal functions (EOFs) and closely related methods (e.g., [4]) find signals that explain maximum variance but do not address the problem of mixing. An appropriate method that addresses the latter problem is independent component analysis (ICA). ICA is a method for separating mixed signals into their sources. It is based on the assumption that different sources, stemming from different physical and dynamical processes, can be considered to be statistically independent and can be separated by ICA.

ICA can be relevant to climate research to identify and separate possible anomaly patterns that may have different forcing [5]. ICA has been used in climate research only lately. For example, [6] applied ICA to the study of the West African vegetation on intraseasonal and interannual time scales. Reference [7] analyzed sea level pressure via ICA to identify the main contributors to the Arctic Oscillation signal. Other authors applied ICA to the sea surface temperature (SST) and surface temperature [810].

Conventional ICA is essentially a one-dimensional algorithm in that the sources are unidimensional. In reality, however, and given the high nonlinear character of climate one cannot expect all source components (i.e., unidimensional) to be statistically independent. Rather, this can be relaxed in favor of finding groups of hidden (independent) sources. This has led to the development of an extension of ICA to the multidimensional ICA also known as Independent Subspace Analysis or ISA [1114], already used in the analysis of the atmospheric variability [15]. Most common algorithms used to obtain independent sources are based on maximizing various information measures such as negentropy (NE).

The SST field is an important source of information to weather and climate variability as it constitutes the bottom boundary forcing of the overlying atmosphere and controls air-sea interaction through sensible and latent heat fluxes. SST variability is characterized by high nonlinearity and non-Gaussianity [1, 1619]. In this manuscript, we analyze the SST using ISA to identify groups of independent sources. We use high-order cumulant tensors of coskewness and cokurtosis combined with high-order singular value decomposition (HOSVD) and minimize the generalized mutual information (MI) to obtain the independent sources. However, the performance of the methods may be hindered by sample sizes, serial correlation, or high dimensionality, well known characteristics of weather/climate. Therefore, we present statistical tests examining the biases of NE and MI and the effects of the sampling and dimensionality, which helps in finding robust non-Gaussian subspaces.

The manuscript is organized as follows. Section 2 presents the separation method and corresponding statistical tests. Section 3 presents application to the SST anomaly data and to a synthetic non-Gaussian surrogate consistent with the same data. A summary and conclusion are provided in the last section.

2. The Method

2.1. Cumulant-Based Independent Subspace Analysis

We designate by our space-time data matrix anomaly (of rank ), representing realizations of a random vector (). Principal Component Analysis (PCA) or Singular Value Decomposition (SVD) [20] are used to construct the PCs: with variances: . The method described below is applied to any whitened or isotropic random vector of dimension obtained from or to a subsequent orthogonal rotation () (note that since is isotropic and standard multivariate Gaussian, then the same holds for ). The source separation, or the factorization of the pdf of , into partial pdfs is not a trivial problem when is non-Gaussian. In fact, this is the case when the joint negentropy (NE), [21], which is invariant under any linear homeomorphism, is strictly positive, that is, where is the multinormal probability density function (pdf) with the same first two moments as .

The pdf separation relies on the maximization of the sum of negentropies of a set of candidate sources. Precisely, let us consider the partition of the rotated vector into a number of subvectors or candidate sources (e.g., scalars, dyads, or tuples) as with . Thanks to the separation theorem [13] or Negentropy Lemma [15], decomposes as a sum of partial negentropies plus the nonnegative generalized mutual information (MI) [22] as The goal of ISA is to find a rotated vector which is split into the largest number of less statistically dependent subvectors minimizing MI. ICA is a restriction of ISA to the case of scalar sources, that is, for all . Gaussian distributed sources, (if they exist) are necessarily unidimensional and span the so-called Gaussian subspace of dimension . Its orthogonal complement (with respect to the covariance inner product), of dimension , is , that is, .

The explicit MI dependence on (with components ) is quite difficult to implement in general. Here, we use a cumulant-based contrast function, relying on the Edgeworth expansion approximation of the pdf of [23], namely: , where is a multivariate linear combination of Hermite polynomials depending on cumulants up to a given maximum order, set here to . The third and fourth-order cumulants of (components of the coskewness and cokurtosis tensors) are for , , , in . They are invariant under permutation of indices (hypersymmetry). After applying a tensor change of basis, the resulting cumulants are written in terms of entries for , , , and in . Cumulant tensors of any order aggregate measures of the joint non-Gaussianity of all possible subvectors of . They vanish for all when is multinormal or for those intersecting independent subvectors. To approximate NE, we introduce a positive defined matrix by contracting and asand similarly, for , that is, . The Edgeworth-based NE approximation is given bywhere, for example, stands for the Frobenius’ squared norm (sum of all squared components) of the tensor . Note that there are examples of non-Gaussian pdfs verifying , including isotropic pdfs, for which nonlinear methods must be used [24, 25], such as the nonlinear-PCA (NL-PCA) [26] which has been applied to climatic data [2731].

The approximation (5) is scaled as , if behaves like an average of an equivalent number of independent and identically distributed (iid) random variables. is appropriate to be used in ISA since it satisfies to an equivalent separation theorem (1):where collects all the coskewness and cokurtosis mixing components of the different candidate sources. Therefore, ISA holds when is minimum over the manifold of orthogonal rotations in . This is related to the Joint Approximated Diagonalization Eigen-matrices’ (JADE) method [32].

The rotation matrix composed of the Singular vectors (SVecs) of (in columns) provides the Tucker decomposition [33] or high-order singular value decomposition (HOSVD) [34]. The HOSVD applied to cumulant tensors lead to the Principal Cumulant Component Analysis (PCCA) [35] and to the Joint Moment Component Analysis (JMCA), [36]. Those techniques find an optimal reduced basis accounting for most of the variance and negentropy explained by a set of cumulants, making thus a generalization of PCA.

The singular values (SVs) of can be expressed as which merge all the negentropy terms where the projection of on the th SVec: takes part. It provides the interpretation for the leading SVecs as being the main axes where non-Gaussianity is sought. Note also that there are other scalar tensor invariants characterizing the joint non-Gaussianity, relying on linear combinations of , (e.g., determinant of ).

The vectors corresponding to nonnull SVs span the estimated non-Gaussian subspace whereas the kernel spanned by the trailing SVecs associated with zero SVs spans the estimated Gaussian subspace as far as third- and fourth-order cumulants are concerned. Note, however, that the method provides an overestimation of the Gaussian manifold and an underestimation of the non-Gaussian manifold .

To summarize, the first step of source separation is to discard the Gaussian subspace and restrict the analysis to vectors in of dimension . Let us denote the non-Gaussian part of by .

A given partition that leads to separated true non-Gaussian sources, that is, , makes the invertible matrix block-diagonal with blocks. The SVD of yields blocks of SVecs: ; , spanning the same source’s subspaces (i.e., ; ) (up to rotation, permutation, or inversion of axes within each source). This, however, does not provide information on how to group the SVecs onto sources (source configuration indeterminacy), which is a complex problem.

In the presence of nonseparable multivariate sources, one will still be interested in finding, within each source, axes maximizing the single NEs. In other words, one solves the ICA problem in by maximizing the NE sum in (6), that is, solving an ICA problem [37, 38] providing independent components (ICs) .

For a large class of pdfs, [13] has showed and conjectured that the ISA solution results in general from the grouping of ICs into independent subspaces.

The general problem of grouping ICs coming from HOSVD followed by ICA benefits from a special decomposition of into positive terms, hereafter called interactivities (ITIs). The decomposition runs through all possible nonnull subsets with dimension , where . Therefore,As an example, if and keeping in mind the hypersymmetry, we get The ITIs, for 1D, that is, (self-interactivities), coincide with single NEs whereas for they coincide with joint-interactivities, which provide a measure of the non-linear statistical relationships within contributing to . Using cumulant properties [39], it is straight forward to show that when has non-null intersections with different independent sources. Therefore, after computing all ITIs, the source identification consists of determining a partition of composed of subvectors of minimal dimension such that for all intersecting at least two subvectors and , . We note that the NE decomposition is not unique. In fact, the concept of interaction information [4042] provides a general frame for the decomposition of in terms of interaction information, defined for every , which are given by linear combinations of mutual information of subvectors . See [43] for a proper application of ITIs on nonlinear fluid dynamics for the assessment of nonlinear wave resonances.

2.2. Independent Subspace Analysis from Finite Samples
2.2.1. Statistical Tests of Multivariate Non-Gaussianity

The above source separation method depends crucially on the approximation of the expectation of any function of -function: (e.g., the cumulants (2), the HOSVD matrix (4), the negentropy (5), and its interactivities (8)), through its estimation obtained by sample averages: . The sample high-order cumulant estimators may be quite biased due to the eventual presence of single and joint outliers on short samples. This is particularly relevant to weather/climate with relatively short datasets with nonnegligible serial correlations yielding a smaller effective number of degrees of freedom: . The bias becomes even larger when we deal with high-order cumulants. In fact, the th order (hypersymmetric) cumulant tensor in dimensions has independent entries, which scales as a power law of . This can lead to a high bias of the cumulants of the leading SVs of and to a possible fictitious non-Gaussianity of the pdf, resulting from the “curse of dimensionality.” Briefly speaking, the “false” non-Gaussianity may be due to either small sample size, high serial correlation, high dimension , or high-order statistics . To avoid non-Gaussianity overestimation, some procedures are usually applied including (i) data projection onto low-dimensional subspaces of high relevance (projection pursuit (PP) rationale [44, 45]); (ii) outlier’s removal; and (iii) use of resistant (or order) statistics considering reasonable values of the order .

Another strategy that we will also follow here is to evaluate the distribution, moments, and quantiles of the distribution of the estimators , necessary to perform ISA, under the null hypothesis () of multivariate Gaussianity of . All quantities under will hereafter be noted with the subscript “.” For iid realizations, all the above statistics are scaled as functions of and , for which one has obtained analytic asymptotic formulas. In the presence of serial correlation, the statistics are derived from a Monte-Carlo procedure by taking an ensemble of surrogate Gaussian time-series with the same , taken here as simple red noise AR processes fitted to each component as ; , where is a standard Gaussian white noise process and is the lag-one autocorrelation. For the significance tests of non-Gaussianity, the upper quantiles (e.g., 90%–99%) statistics are taken from a large ensemble of surrogates.

The first step of source decomposition is to test the null hypothesis . The finite-sample estimator of , hereafter denoted as (all quantities based on sample statistics are overlined with a tilde) and others based on truncated NE approximations are conservative tests of (or less powerful tests of its negation: ~), because they reject non-Gaussianity (Type II error) resulting from discarded cumulants of order .

Now, we discuss the behavior of the estimator of (under ) computed from -sized samples of a vector of independent scalars. The pdf of is positively skewed, depending on the joint distribution of and , which are in general not independent since they satisfy a number of inequality constraints (e.g., ). The bias of is given bywhere the first and second square-bracket terms of (10) represent, respectively, scaled versions of and . The bias (10) tends asymptotically to . It also provides a lower-bound of the positive bias for non-Gaussian pdfs, as it occurs for any negentropy estimator [46].

The biases from the Frobenius norms were obtained by multiplying the biases of the different squared cumulant estimators by their multiplicities [36]. For instance, there are cumulant biases of type ; in based on the moments of the standard Gaussian: . Figure 1 shows the dependence on and for , of the ratiobetween the actual bias and the asymptotic bias of (at ), obtained from 1000 Monte-Carlo replicates. Moreover, the dependence of on is quite negligible with relative deviations less than 5% (not shown).

From the analysis of Figure 1(a) and as expected, the ratio remains close to 1 when is small and is large. For , the ratio has reached the asymptotic value, where it becomes dependent only on with , where is interpreted as the “average interval between independent realizations.” For small sample sizes, , due to a negative higher-order term in (10).

The variance of involves covariances between generic cumulants, for example, and , of subvectors . Those covariances satisfy ; with if and if and otherwise.

As for the bias, the quantiles of hereby noted scale asymptotically as . Figure 1(b) shows the dependence on and for , of the ratio . This figure provides a kind of a rule of thumb thresholds of for rejecting in the asymptotic regime.

2.2.2. Estimation of the Gaussian Subspace by HOSVD

Following a rejection of , the next step is to estimate the Gaussian subspace and its dimension . With finite samples, one computes the sample-based matrix of (4). However, can no longer be estimated by (see Section 2.1) but through the subspace spanned by SVecs associated with the trailing (not significant) SVs.

Under , the (estimated) SVs , (see (7)), with quantiles , , are asymptotically and the SVs become undistinguishable for large [36]. However, the SVs are not independent and consequently their marginal quantiles cannot be estimated independently for each . To overcome this difficulty, a sequence of statistical tests is applied, as discussed in the next subsection.

In summary, we first consider the vectors and , where is a generic surrogate of and and are matrices filled by the ranked SVecs of and , respectively. For each in , we construct the vector , where the trailing components are swapped by the associated ranked surrogates. Then we compute their NE contributions (5) and finally the quantile evaluated from a large ensemble of surrogates. The test then proceeds backward, starting from the smallest SV (). Therefore, if , for all , then the null hypothesis (undistinguishable from Gaussian surrogates for all ) cannot be rejected. The smallest for which is not rejected determines the estimated Gaussian subspace with dimension (spanned by the trailing SVecs of ), and (spanned by the leading SVecs of ).

2.2.3. Estimation of Multivariate Sources

The next step is to identify the non-Gaussian sources within the vector of non-Gaussian components spanning . Since is a subvector of , we easily get .

We test the separability of into scalar and/or subvector sources by looking for a new orthogonal rotation composed of a maximal number of subvector candidate sources , , such that the joint pdf satisfieswhere the norm of the pdf separation error is below the threshold of rejection of the null hypothesis of nonseparability (). There is an a priori source configuration indeterminacy. For instance, for , there are five possible ways, namely: 4 scalars (ICA case); 2 scalars and a dyad; 2 dyads; a triad plus a scalar and the global quartet.

For a chosen source configuration, an appropriate contrast function to maximize iswhere must be below some significance threshold, passing the separation statistical test. However, an undecidable situation may occur when multiple configurations pass the test. To avoid the case of multiple configuration indeterminacy, we solve instead the ICA problem providing the vector of best scalar sources: . The algorithm to get is explained in Appendix A. Finally, we build multivariate sources forming appropriate groups of components, according to the conjecture of [13], and as below.

We start with the NE decomposition into interactivities (8) asTo compare and account for the statistical significance of -interactivities in (15) with different dimensions , we introduce scaled interactivities through its comparison with the quantiles of the NE of the ranked (by HOSVD) Gaussian surrogates of . It is is defined asWhen in (16), is the rotation of towards isotropy (e.g., by Gram Schmidt orthogonalization) with being composed of the standard Gaussianized components of . Since has Gaussian marginals, this procedure leads to a fairer and unbiased test when has non-Gaussian pdf marginals. This is also supported by the fact that the mutual information of is where the second term is the Gaussian MI component [47, 48], which depends on the correlation matrix of and the third term is the sum of the NEs of the components of . The above measure is therefore resistant to single outliers of the marginal pdfs.

The different scaled interactivities (16) are then sorted by decreasing order, retaining uniquely the leading ones verifying where is a fixed level of acceptance. The emerging large interactivities normally yield ICs, which contribute to the multivariate sources. The NE interactivities are then grouped along the relevant interactivities until the stopping criteria are reached, for example, , which will be adopted in practice.

3. Analysis of Sea Surface Temperature Data

3.1. Data Description and Processing

We consider the monthly SST data from the “Extended Reconstruction Sea Surface Temperature” (ERSST) version v3b product (, from the International Comprehensive Ocean–Atmosphere Data Set (ICOADS) release 2.4. The SST data are on latitude-longitude resolution and span the period 1875 to present [49]. We only consider the period 1910–2011 (102 years), restricted to the region equatorward of 65°. At each grid point, the climatology and the (12-month moving average) linear trend in addition to the mean seasonal cycle are removed from the data [50], thus yielding monthly SST anomalies (SSTAs), on grid points. A EOF analysis is then applied [51] and we only use the space spanned by the leading mode for the application of ICA/ISA. Its output is independent from PC variances, even when the used PCs are quasi-degenerated. The total variance is (°C)2, the local variance at point is , with an average ~(0.50)2 (°C)2 ( = average of ) and we let the fractional local variance . Figure 2(a) and Table 2 show the % of explained variance of the 12 leading EOFs with their confidence intervals [52] and Figure 2(b) shows the fitted lag-one autocorrelation based on an AR model.

Figure 3 shows the leading 12 EOFs which are described, with references in Table 1. Based on , an estimate of the “average interval between independent realizations” is , (Figure 1(a)) which is not too different from (e.g., [53]).

EOF rankDescription with the main SSTA modesAvailable index PC/Index cor.Refs.Notes

1El Niño (ENSO)MEI0.98[56]

2Atlantic Multi-Decadal Oscillation (AMO)AMO0.80[57]

3Pacific Decadal Oscillation (PDO): negative phasePDO[58, 59]

4Northern Pacific Gyre Oscillation (NPGO)NPGO[60]
El Niño Modoki (EM): negative phaseEMI[61]Also known as Central Pacific or noncanonical El Niño

5El Niño Modoki (EM): positive phaseEMI[61]
Subtropical Indian Ocean Dipole: negative phaseSIOD[62, 63]First EOF mode of the SSTA in the Indian Ocean ranging in the 15°S–45°S band

6Atlantic Niño: negative phaseATL3[64, 65]
South Atlantic Subtropical DipoleSASD[66]First EOF coupled mode of the SSTA/wind in the South Atlantic
Subtropical Indian Ocean Dipole: positive phaseSIOD[62, 63]First EOF mode of the SSTA in the Indian Ocean ranging in the 15°S–45°S band
Zonally oriented dipole in South Pacific[67, 68]Second EOF of the South Pacific SSTA

7Meridional sequence of dipoles crossing the North-South Pacific
South Pacific Ocean Dipole (SPOD)SPOD[67, 69]First EOF of the South Pacific SSTA

8North Atlantic tripole (NAT)NAT[70]
Tropical Atlantic SST dipole[68]

9South Pacific Ocean Dipole (SPOD)SPOD[67, 69]First EOF of the South Pacific SSTA

10Zonally oriented dipole in South Atlantic[68]
Zonally oriented dipole in South Pacific[67, 68]Second EOF of the South Pacific SSTA

11Meridional Pacific wavelike pattern
Meridional Atlantic wavelike pattern
Subtropical Indian Ocean Dipole (positive phase)SIOD[62, 63]First EOF mode of the SSTA in the Indian Ocean ranging in the 15°S–45°S band

12Tripolar pattern in North Pacific
North Atlantic tripole (NAT) NAT[70]



Figure 2(b) shows that varies from about 3 months (PC10) to more than 12 months (PC2). For our simulation later, we will use the geometric average of the 12 values, that is, 0.86 (see Section 3.3).

3.2. Negentropy Distribution and Non-Gaussian Subspace Estimation

Non-Gaussianity and negentropy of a field have counterparts both on a local and on a PC basis. In fact, local cumulants are obtained from cumulants of the whitened PCs composing vector , through linear operators. Let us consider, for example, the local-basis coskewness tensor . Its Frobenius norm is the area integral of the local squared skewness. By using the EOF orthogonality property, we may express this as a weighted sum over all grid point triads :where stand for the relative area at grid point and is the fractional local variance. The contributions to (17) of a triad of whitened PCs are proportional to their variance fractions. The Frobenius norm (17) is dominated by contributions primarily from regions of high variance and absolute skewness, namely, the positively skewed El Niño zone [19, 54], in addition to to the positive skewness southwards of 50°S [17].

Now, the NE can be decomposed using (5). For each , we denote by the set of leading whitened PCs. Figure 4 shows the diagonal terms (), Figure 4(a) and the total NE , Figure 4(b). Note that, for a given , the smaller values of are either the least non-Gaussian PCs or the less interactive PCs with the remaining PCs (), providing therefore good candidates for the Gaussian subspace of . For example, with , we have , and 12 (Figure 4(a)). Statistical significance of NE components and totals are assessed by using quantiles computed from Gaussian surrogates (Figure 4). Notice, for instance, the confidence level <70% of of the PC2 due to its small degrees of freedom.

Biases may be approximated as (see (11)) with , leading to (see Figure 1(a)). For instance, for , , the bias is 1.7 (55%) with the remaining value (1.4) being attributed to true non-Gaussianity. Therefore, one expects that biases shall be drastically reduced after discarding the Gaussian subspace.

Figure 4(a) shows that, for all , the largest NE component is , resulting mostly from the nonlinear correlations of PC1 with other PCs and from its self-NE, due to its positive skewness (Table 2).

Using the nonnegative consecutive jumps: (Figure 4(a)), one may detect which interactivities mostly contribute to . The largest differences are for pairs: , in consistency with the high confidence level of for . Those jumps are often due to particular joint cumulants in the 2D subspaces spanned by , responsible for particular types of joint non-Gaussianity, for example, scatter-plots resulting for instance from the occurrence of joint extremes. For instance, , is favored for the sign combinations and of . As an illustration, Figures 5(a) and 5(b) show, respectively, scatter-plots of and (the marginals have been Gaussianized (subscript g) to emphasize the 2D pdf difference with respect to 2D isotropic Gaussianity). Figure 5 shows that most of the non-Gaussianity is due to the occurrence of joint coherent 2D positive extremes, for example, the 1983 and 1997 extreme El Niño events (see Figure of [50]) during which a positive PC11 and a negative PDO phase have occurred.

Since is quite low and not highly significant (Figure 4(a)), we have decided to apply ICA and ISA in the 11D subspace (see Section 3.4) in which the total NE reach a confidence level > 95%. Besides the test to limit [55] presents an alternative limit. The dimensionality reduction, however, is performed below via the statistical detection of the non-Gaussian subspace as detailed in Section 2.2.2.

Therefore, we now consider the state vector and the associated vector , yielding , and a bias of 1.73, leaving ~0.82 for “true non-Gaussianity.” Figure 6 shows the NE components (or SVs) (), along with the threshold quantiles 90% and 95% of obtained from an ensemble of 1000 Gaussian surrogates. For comparison, we also show the SVs (open squares) of surrogate realizations of . The SVs of all lie below the above thresholds, supporting the consistency of the Gaussianity test. Figure 6 also suggests that the estimated Gaussian subspace of the SSTA is spanned by the last 6 SVecs (with a confidence level ) and so and, for the non-Gaussian subspace . The total NE restricted to is with a bias of 0.14. The “useful” true non-Gaussianity is 0.97 which is of the order of the above estimated value (0.82). The NE components (open black circles), with their low values compared to , , are also shown in Figure 6.

The squared projections of onto each of the 11 whitened PCs are given by the inner product (Appendix B) , for , (Table 2), where is the (unitary) th vector. By definition, those values add up to . Table 2 shows that some PCs mostly project onto the non-Gaussian subspace such as PC1 (El Niño), PC4 (NPGO), and PC10 while others mostly project onto the Gaussian subspace, for example, PC7, PC8, and PC3 (negative PDO). The 48% of explained variance is distributed, respectively, as = 24.7% and 23.3%, for the non-Gaussian and Gaussian subspaces.

3.3. Application to Synthetic Non-Gaussian Data

We generate AR independent series , , . The whitened and uncorrelated components of the non-Gaussian vector are computed followingwith measuring the “strength” of non-Gaussianity. So is composed of Gaussian scalar sources spanning and non-Gaussian variables spanning . The latter is composed of one dyad and a non-Gaussian triad .

The next step consists of testing the null hypothesis of global Gaussianity (i.e., ) and evaluating by varying in and in .

Figure 7(a) shows the Monte-Carlo average obtained from 1000 surrogates. As expected, the mean NE increases with (nonlinearity effect) and (dimensionality effect), though at a much smaller rate. The bias is the value reached at in agreement with (11).

3.3.1. Gaussian Subspace Recovery

We investigate the strength of as a function of and . The null hypothesis is rejected at the confidence level if using (Figure 1(b)). The Monte-Carlo fraction of rejection is shown in Figure 7(b) for . The contour line at provides the graph of for the recovery of the non-Gaussian sources. As expected, increases with , starting at and at and reaching at . This means that the larger , the stronger (through ) the non-Gaussian embedded sources. Moreover, diminishes under a less conservative test (lower ) or when the effective number of temporal degrees of freedom increases. This puts a limit to what is recoverable from the available (often short) weather/climate data, and, therefore, less conservative tests of Gaussianity or data dimension reduction must be considered (e.g., ) [47].

3.3.2. ICA and ISA

Here, we apply ICA and ISA for () and compute the skill of source separation in two cases: (a) , high noise, and (b) , intermediate noise, yielding, respectively, total NE of 1.79 and 2.57.

Figures 8(a) and 8(b) are similar to Figure 6 but for the cases (a) and (b), respectively. The analysis yields the correct dimension in both cases: with despite the fact that the SVs in the “high noise” case (a) are smaller. The score of the non-Gaussian subspace recovery is assessed using the inner products (compared to 5 in the ideal recovery). Those products are, respectively, 3.99 and 4.49 for cases (a) and (b). Lower values of lead to a subestimation of and to a worse estimation of .

Now, by applying ICA (see Section 2.2.3 and Appendix A) we get ICs making the vector . Tables 3 and 4 give, respectively, for cases (a) and (b), a summary of ICA with the sorted maximized self-NEs , the NE components , and the IC vector loadings in . A few conclusions can be drawn from these tables; (i) the sum of self-NEs of ICs (2nd column) cannot explain the total NE (3rd column). This corroborates the existence of NE originating from IC interactivities and (ii) Up to a certain statistical significance, a given IC projects uniquely onto a certain source as seen through the dominant projections of ICs (in bold font). From these tables, the estimated dyad is spanned by IC1 and IC4 whereas the estimated triad is spanned by IC2, IC3, and IC5. That supports the conjecture of [13] about the source emerging by IC grouping. Recovered sources may be rotated or inverted with respect to initial data. The distribution of ICs among sources may differ for smaller . The squared projections onto the true sources are, respectively, 1.35 (on a maximum of 2) and 2.42 (on a maximum of 3) in case (a) and 1.68 and 2.75 in case (b). Like the non-Gaussian subspace, the recovery of individual sources is more skillful in case (b).

IC rank (dyad) (dyad) (triad) (triad) (triad)

IC1: 0.1270.162−
IC2: 0.0660.125−0.070.07
IC3: 0.0420.102−0.01−0.08
IC4: 0.0140.0530.210.290.01
IC5: 0.0110.063−0.030.02


IC rank (dyad) (dyad) (triad) (triad)