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

Authors: Carlos A. L. Pires; Abdel Hannachi
Journal: Complexity
 (2017)

---

## Abstract

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.


---

## Body

## 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 [1–3]. 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 [8–10].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 [11–14], 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, 16–19]. 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. Section2 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 byX∈RNp×Nt our (Np×Nt) space-time data matrix anomaly (of rank Nr), representing Nt realizations of a random vector x(t)∈RNp (t=1,…,Nt). Principal Component Analysis (PCA) or Singular Value Decomposition (SVD) [20] are used to construct the PCs: xPC(t)≡XPC,1t,…,XPC,Nrt′ with variances: λ1≥⋯≥λNr. The method described below is applied to any whitened or isotropic random vector y of dimension N≤Nr obtained from xPC(t) or to a subsequent orthogonal rotation z=R′y (RR′=1N) (note that since y is isotropic and standard multivariate Gaussian, then the same holds for z). The source separation, or the factorization of the pdf of y, ρy into partial pdfs is not a trivial problem when y 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, J(z)≡E[log(ρz/φz)]=J(y)≡Jrot>0 where φz is the multinormal probability density function (pdf) with the same first two moments as z.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 vectorz into a number Ns of subvectors or candidate sources (e.g., scalars, dyads, or tuples) as z′=[z1′,z2′,…,zNs′]=(Z1,…,ZN) with ∑i=1Nsdim⁡zi=N. Thanks to the separation theorem [13] or Negentropy Lemma [15], Jrot decomposes as a sum of partial Ns negentropies plus the nonnegative generalized mutual information (MI) I(z(1),z(2),…,z(Ns)) [22] as (1)Jrot⁡-∑i=1NsJzi=Iz1,z2,…,zNs≡Elog⁡ρz∏i=1Nsρzi≥0,Jzi=12dim⁡zilog⁡2πe+Elog⁡ρzi≥0.The goal of ISA is to find a rotated vector z=R′y which is split into the largest number Ns of less statistically dependent subvectors minimizing MI. ICA is a restriction of ISA to the case of scalar sources, that is, dim⁡(z(i))=1 for all i. Gaussian distributed sources, (if they exist) are necessarily unidimensional and span the so-called Gaussian subspace Wg(z) of dimension Ng. Its orthogonal complement (with respect to the covariance inner product), of dimension Nng=N-Ng, is Wng(z)=Wg(z)⊥, that is, Wg(z)⊕Wng(z)=RN.The explicit MI dependence onR (with components Rij,…,i,j=1,…,N) 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 y [23], namely: ρy≈φ(y)[1+PEd(y)], where PEd is a multivariate linear combination of Hermite polynomials depending on cumulants up to a given p=pmax maximum order, set here to pmax=4. The third and fourth-order cumulants of y (components of the coskewness and cokurtosis tensors) are (2)Syijk≡EYiYjYk;Kyijkl≡EYiYjYkYl-EYiYjEYkYl-EYiYkEYjYl-EYiYlEYjYkfor i, j, k, l in {1,…,N}. They are invariant under permutation of indices (hypersymmetry). After applying a tensor change of basis, the resulting cumulants are (3)Szijk≡∑p,q,r=1NRpiRqjRrkSypqr;Kzijkl≡∑p,q,r,s=1NRpiRqjRrkRslKypqrs,written in terms of R entries for i, j, k, and l in {1,…,N}. Cumulant tensors of any order p≥3 aggregate measures of the joint non-Gaussianity of all possible subvectors w of y. They vanish for all w when y is multinormal or for those w intersecting independent subvectors. To approximate NE, we introduce a positive defined matrix by contracting S and K as(4)Myij≡112∑p,q=1NSyipqSyjpq+148∑p,q,r=1NKyipqrKyjpqr,and similarly, for z, that is, M(z)ij=∑p,q=1NRpiRqjM(y)pq. The Edgeworth-based NE approximation is given by(5)JEdy≡TrMy=JEdz=∑i=1NMiiz=112SzF2+148KzF2,where, for example, S(z)F2 stands for the Frobenius’ squared norm (sum of all squared components) of the tensor S. Note that there are examples of non-Gaussian pdfs verifying JEd(y)=0, 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 [27–31].The approximation (5) is scaled as Jrot-JEd(z)=O(Neq-1), if z behaves like an average of an equivalent number Neq of independent and identically distributed (iid) random variables. JEd(z) is appropriate to be used in ISA since it satisfies to an equivalent separation theorem (1):(6)JEdz=∑i=1NsJEdzi+IEdz1,z2,…,zNs,where IEd collects all the coskewness and cokurtosis mixing components of the different candidate sources. Therefore, ISA holds when IEd is minimum over the manifold SO(N) of orthogonal rotations in RN. This is related to the Joint Approximated Diagonalization Eigen-matrices’ (JADE) method [32].The rotation matrixR=(r1,r2,…,rN) composed of the Singular vectors (SVecs) of M(y) (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) ofM(y):η1,N≥η2,N≥⋯≥ηN,N can be expressed as (7)ηi,N=112∑p,q=1NSzipq2+148∑p,q,r=1NKzipqr2;1≤i≤N,which merge all the negentropy terms where the projection of z on the ith SVec: Zi=ri′y 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 Tr[Mk(z)], k≤N (e.g., determinant of M).The vectors corresponding to nonnull SVs span the estimated non-Gaussian subspaceWng-Ed(z) whereas the kernel Ker(M) spanned by the trailing SVecs associated with zero SVs spans the estimated Gaussian subspace Wg-Ed(z) as far as third- and fourth-order cumulants are concerned. Note, however, that the method provides an overestimation of the Gaussian manifold Wg-Ed(z)⊇Wg(z) and an underestimation of the non-Gaussian manifold Wng-Ed(z)⊆Wng(z).To summarize, the first step of source separation is to discard the Gaussian subspace and restrict the analysis to vectors inWng-Ed(z)=Ker[M(z)]⊥ of dimension Nng-Ed. Let us denote the non-Gaussian part of z by zng=Z1,…,ZNng-Ed′.A given partitionzng′=[z1′,z2′,…,zNs′] that leads to separated true non-Gaussian sources, that is, IEd(z(1),z(2),…,z(Ns))=0, makes the invertible matrix M(zng) block-diagonal with Ns blocks. The SVD of M(zng) yields blocks of SVecs: w(i); i=1,…,Ns, spanning the same source’s subspaces (i.e., Span[w(i)]=Span[z(i)]; i=1,…,Ns) (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 inWng-Ed by maximizing the NE sum ∑i=1Nng-EdJEd(Zi) in (6), that is, solving an ICA problem [37, 38] providing independent components (ICs) Zi.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 ofJEd(zng) into positive terms, hereafter called interactivities (ITIs). The decomposition runs through all possible nonnull Nint≡2m-2m-min(m,pmax) subsets u⊆zng with dimension 1≤|u|≤min(m,pmax=4), where (m≡Nng-Ed). Therefore,(8)JEdzng≡∑u⊆zJEd-ITu,JEd-ITu≡112∑i,j,k:ui,uj,uk=uSuijk2+148∑i,j,k,l:ui,uj,uk,ul=uKuijkl2.As an example, if u=U1,U2′ and keeping in mind the hypersymmetry, we get (9)JEd-ITu=1123S1122+3S2212+1484K11122+6K11222+4K22212.The ITIs, for 1D, that is, dim(u)=1 (self-interactivities), coincide with single NEs whereas for dim(u)≥2 they coincide with joint-interactivities, which provide a measure of the non-linear statistical relationships within u contributing to JEd(zng). Using cumulant properties [39], it is straight forward to show that JEd-IT(u)=0 when u has non-null intersections with different independent sources. Therefore, after computing all ITIs, the source identification consists of determining a partition {u(1),u(2),…,u(Ns)} of z composed of subvectors of minimal dimension such that JEd-IT(v)=0 for all v intersecting at least two subvectors u(α) and u(β), α≠β. We note that the NE decomposition is not unique. In fact, the concept of interaction information [40–42] provides a general frame for the decomposition of J(zng) in terms of interaction information, defined for every u⊆zng, which are given by linear combinations of mutual information of subvectors v⊆u. 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 functionf(z) of z-function: E[f(z)] (e.g., the cumulants (2), the HOSVD matrix M(z) (4), the negentropy (5), and its interactivities (8)), through its estimation obtained by sample averages: f(z)¯. 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: Ndof≤Nt. The bias becomes even larger when we deal with high-order cumulants. In fact, the pth order (hypersymmetric) cumulant tensor in N dimensions has CpN+p-1=(N+p-1)!/[(N-1)!p!] independent entries, which scales as a power law of N. This can lead to a high bias of the cumulants of the leading SVs of M(z) 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 N, or high-order statistics p. 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 p.Another strategy that we will also follow here is to evaluate the distribution, moments, and quantiles of the distribution of the estimatorsf(z)¯G, necessary to perform ISA, under the null hypothesis (HG) of multivariate Gaussianity of z. All quantities under HG will hereafter be noted with the subscript “G.” For iid realizations, all the above statistics are scaled as functions of Nt and N, 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 Ndof, taken here as simple red noise AR(1) processes fitted to each z component as g(t+1)=c^1g(t)+(1-c^12)1/2w(t); t=0,…,Nt, where w(t) is a standard Gaussian white noise process and c^1 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 hypothesisHG. The finite-sample estimator of JEd, hereafter denoted as J~Ed (all quantities based on sample statistics are overlined with a tilde) and others based on truncated NE approximations are conservative tests of HG (or less powerful tests of its negation: ~HG), because they reject non-Gaussianity (Type II error) resulting from discarded cumulants of order p>pmax=4.Now, we discuss the behavior of the estimatorJ~Ed:G of JEd:G (under HG) computed from Nt-sized samples of a vector of N independent scalars. The pdf of J~Ed:G is positively skewed, depending on the joint distribution of S~G and K~G, which are in general not independent since they satisfy a number of inequality constraints (e.g., K1111-S1112+2≥0). The bias of J~Ed:G is given by(10)EJ~Ed:G=1NdofPJEdN+O1Ndof2,1NdofPJEdN≡1121NdofN3+3N2+2N+1481NdofN4+6N3+11N2+6N,where the first and second square-bracket terms of (10) represent, respectively, scaled versions of E[S~GF2] and E[K~GF2]. The bias (10) tends asymptotically to PJEd(N)/Ndof. It also provides a lower-bound of the positive J~Ed 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 3N(N-1) cumulant biases of type E[(YiYiYj)¯2]=2Ndof-1+O(Ndof-2); i≠j in E[S¯GF2] based on the moments of the standard Gaussian: E(U2k)=(2k)!/(2kk!). Figure 1 shows the dependence on Nt and c^1 for N=4, of the ratio(11)fREdNt,c^1,N≡EJ~Ed:GNtPJEdNbetween the actual bias and the asymptotic bias of J~Ed:G (at c^1=0), obtained from 1000 Monte-Carlo replicates. Moreover, the dependence of fREd on N is quite negligible with relative deviations less than 5% (not shown).Figure 1
(a) Contour plot offREd as a function of Nt, c^1 for a fixed dimension N=4. (b) Ratio fREd-Psig for the confidence level psig=0.95, as a function of the dimension N and the lag-one autocorrelation c^1 for a sample size Nt=1000. The ratio increases with decreasing N or increasing c^1.
(a)
(b)From the analysis of Figure1(a) and as expected, the ratio fREd remains close to 1 when c^1 is small and Nt is large. For Nt≥500, the ratio has reached the asymptotic value, where it becomes dependent only on c^1 with Ndof=Nt/fREd(∞,c^1,N)<Nt, where fREd is interpreted as the “average interval between independent realizations.” For small sample sizes, fREd<fREd(∞,c^1,N), due to a negative higher-order term in (10).The variance ofJ~Ed:G involves covariances between generic cumulants, for example, cm(u1) and cm(u2), of subvectors u1,u2⊆z. Those covariances satisfy cov[cm(u1),cm(u2)]=kNt-d; k≥0 with d=2 if u1=u2 and d=3 if u1⊃u1∩u2≠⌀ and d=∞ otherwise.As for the bias, the quantilespsig∈[0,1] of J~Ed:G hereby noted Qpsig(J~Ed:G) scale asymptotically as Nt-1. Figure 1(b) shows the dependence on N and c^1 for Nt=1000, of the ratio fREd-Psig≡Qpsig(J~Ed:G)[Nt/PJEd(N)]. This figure provides a kind of a rule of thumb thresholds of J~Ed for rejecting HG in the asymptotic regime.


#### 2.2.2. Estimation of the Gaussian Subspace by HOSVD

Following a rejection ofHG, the next step is to estimate the Gaussian subspace Wg and its dimension Ng. With finite samples, one computes the sample-based matrix M~(z) of M(z) (4). However, Wg can no longer be estimated by Ker(M~) (see Section 2.1) but through the subspace spanned by SVecs associated with the trailing (not significant) SVs.UnderHG, the (estimated) SVs η~i,N:G, (see (7)), with quantiles Qpsig(η~i,N:G), i=1,…,N, are asymptotically O(Ndof-1) and the SVs become undistinguishable for large Ndof [36]. However, the SVs are not independent and consequently their marginal quantiles cannot be estimated independently for each i. To overcome this difficulty, a sequence of statistical tests is applied, as discussed in the next subsection.In summary, we first consider the vectorsz=Z1,…,ZN′=R′y and zs≡Zs,1,…,Zs,N′=Rs′ys, where ys is a generic surrogate of y and R and Rs are matrices filled by the ranked SVecs of M~(y) and M~(ys), respectively. For each i in {1,…,N}, we construct the vector bi≡(Z1,…,Zi-1,Zs,i,…,Zs,N), where the trailing i,…,N components are swapped by the associated ranked surrogates. Then we compute their NE contributions M~ii(bi) (5) and finally the quantile Qpsig[M~ii(bi)] evaluated from a large ensemble of surrogates. The test then proceeds backward, starting from the smallest SV (i=N). Therefore, if M~kk(z)<Qpsig[M~kk(bk)], for all k≥i, then the null hypothesis HG,i (Zkundistinguishable from Gaussian surrogates Zs,k for all k≥i) cannot be rejected. The smallest i for which HG,i is not rejected determines the estimated Gaussian subspace Wg-Ed with dimension Ng-Ed (spanned by the trailing Ng-Ed SVecs of R), and Wng-Ed≡Wg-Ed⊥ (spanned by the leading Nng-Ed≡N-Ng-Ed SVecs of R).


#### 2.2.3. Estimation of Multivariate Sources

The next step is to identify the non-Gaussian sources within the vector of non-Gaussian componentszng≡Z1,…,ZNng-Ed′ spanning Wng-Ed. Since zng is a subvector of z, we easily get J~Ed(zng)≤J~Ed(z).We test the separability ofWng-Ed into scalar and/or subvector sources by looking for a new orthogonal rotation zISA≡RISA′zng=Z1,ISA,…,ZNng-Ed,ISA′=z1,ISA′,z2,ISA′,…,zNs,ISA′′ composed of a maximal number Ns≥2 of subvector candidate sources z(i),ISA, i=1,…,Ns, such that the joint pdf satisfies(12)ρzISA=1+εISA∏k=1Nsρzk,ISA,where the norm |εISA| of the pdf separation error εISA is below the threshold of rejection of the null hypothesis of nonseparability (Ns=1). There is an a priori source configuration indeterminacy. For instance, for Nng-Ed=4, 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 is(13)FISAR≡∑i=1NsJ~Edzi,ISA=J~Edzng-I~Edz1,ISA,z2,ISA,…,zNs,ISA,where I~Ed 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 (14)RICA≡argmaxRISA∑i=1NgJ~EdZi,ISA,providing the vector of best scalar sources: zICA=RICA′zng=Z1,ICA,…,ZNNg-Ed,ICA′. The algorithm to get RICA is explained in Appendix A. Finally, we build multivariate sources forming appropriate groups of zICA components, according to the conjecture of [13], and as below.We start with the NE decomposition intoNint interactivities (8) as(15)J~Edzng=J~EdzICA=∑u⊆zICAJ~Ed-ITu.To compare and account for the statistical significance of u-interactivities in (15) with different dimensions dim⁡(u), we introduce scaled interactivities J~Ed-ITsig(u) through its comparison with the quantiles of the NE J~Ed-IT(uG) of the ranked (by HOSVD) Gaussian surrogates uG of u. It is is defined as(16)J~Ed-ITsigu=J~Ed-ITwQpsigJ~Ed-ITuG;w=uifdimu=1uGsiifdimu>1.When dim(u)>1 in (16), w is the rotation of uGs towards isotropy (e.g., by Gram Schmidt orthogonalization) with uGs being composed of the standard Gaussianized components of u. Since uG has Gaussian marginals, this procedure leads to a fairer and unbiased test when u has non-Gaussian pdf marginals. This is also supported by the fact that the mutual information of u is I(u)=I(uGs)=I(uGsi)+[-1/2log⁡(det⁡CGs)]+∑k=1NJ(uGsi,k)≈I(uGsi) where the second term is the Gaussian MI component [47, 48], which depends on the correlation matrix CGs of uGs and the third term is the sum of the NEs of the components of uGsi. The above measure is therefore resistant to single outliers of the marginal pdfs.The differentNint scaled interactivities (16) are then sorted by decreasing order, retaining uniquely the leading ones verifying J~Ed-ITsig(u)≥Jsig where Jsig 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, ∑uJ~Ed-ITsig(u)<J~Ed(zng)-E[J~Ed(zng)], which will be adopted in practice.


## 2.1. Cumulant-Based Independent Subspace Analysis

We designate byX∈RNp×Nt our (Np×Nt) space-time data matrix anomaly (of rank Nr), representing Nt realizations of a random vector x(t)∈RNp (t=1,…,Nt). Principal Component Analysis (PCA) or Singular Value Decomposition (SVD) [20] are used to construct the PCs: xPC(t)≡XPC,1t,…,XPC,Nrt′ with variances: λ1≥⋯≥λNr. The method described below is applied to any whitened or isotropic random vector y of dimension N≤Nr obtained from xPC(t) or to a subsequent orthogonal rotation z=R′y (RR′=1N) (note that since y is isotropic and standard multivariate Gaussian, then the same holds for z). The source separation, or the factorization of the pdf of y, ρy into partial pdfs is not a trivial problem when y 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, J(z)≡E[log(ρz/φz)]=J(y)≡Jrot>0 where φz is the multinormal probability density function (pdf) with the same first two moments as z.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 vectorz into a number Ns of subvectors or candidate sources (e.g., scalars, dyads, or tuples) as z′=[z1′,z2′,…,zNs′]=(Z1,…,ZN) with ∑i=1Nsdim⁡zi=N. Thanks to the separation theorem [13] or Negentropy Lemma [15], Jrot decomposes as a sum of partial Ns negentropies plus the nonnegative generalized mutual information (MI) I(z(1),z(2),…,z(Ns)) [22] as (1)Jrot⁡-∑i=1NsJzi=Iz1,z2,…,zNs≡Elog⁡ρz∏i=1Nsρzi≥0,Jzi=12dim⁡zilog⁡2πe+Elog⁡ρzi≥0.The goal of ISA is to find a rotated vector z=R′y which is split into the largest number Ns of less statistically dependent subvectors minimizing MI. ICA is a restriction of ISA to the case of scalar sources, that is, dim⁡(z(i))=1 for all i. Gaussian distributed sources, (if they exist) are necessarily unidimensional and span the so-called Gaussian subspace Wg(z) of dimension Ng. Its orthogonal complement (with respect to the covariance inner product), of dimension Nng=N-Ng, is Wng(z)=Wg(z)⊥, that is, Wg(z)⊕Wng(z)=RN.The explicit MI dependence onR (with components Rij,…,i,j=1,…,N) 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 y [23], namely: ρy≈φ(y)[1+PEd(y)], where PEd is a multivariate linear combination of Hermite polynomials depending on cumulants up to a given p=pmax maximum order, set here to pmax=4. The third and fourth-order cumulants of y (components of the coskewness and cokurtosis tensors) are (2)Syijk≡EYiYjYk;Kyijkl≡EYiYjYkYl-EYiYjEYkYl-EYiYkEYjYl-EYiYlEYjYkfor i, j, k, l in {1,…,N}. They are invariant under permutation of indices (hypersymmetry). After applying a tensor change of basis, the resulting cumulants are (3)Szijk≡∑p,q,r=1NRpiRqjRrkSypqr;Kzijkl≡∑p,q,r,s=1NRpiRqjRrkRslKypqrs,written in terms of R entries for i, j, k, and l in {1,…,N}. Cumulant tensors of any order p≥3 aggregate measures of the joint non-Gaussianity of all possible subvectors w of y. They vanish for all w when y is multinormal or for those w intersecting independent subvectors. To approximate NE, we introduce a positive defined matrix by contracting S and K as(4)Myij≡112∑p,q=1NSyipqSyjpq+148∑p,q,r=1NKyipqrKyjpqr,and similarly, for z, that is, M(z)ij=∑p,q=1NRpiRqjM(y)pq. The Edgeworth-based NE approximation is given by(5)JEdy≡TrMy=JEdz=∑i=1NMiiz=112SzF2+148KzF2,where, for example, S(z)F2 stands for the Frobenius’ squared norm (sum of all squared components) of the tensor S. Note that there are examples of non-Gaussian pdfs verifying JEd(y)=0, 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 [27–31].The approximation (5) is scaled as Jrot-JEd(z)=O(Neq-1), if z behaves like an average of an equivalent number Neq of independent and identically distributed (iid) random variables. JEd(z) is appropriate to be used in ISA since it satisfies to an equivalent separation theorem (1):(6)JEdz=∑i=1NsJEdzi+IEdz1,z2,…,zNs,where IEd collects all the coskewness and cokurtosis mixing components of the different candidate sources. Therefore, ISA holds when IEd is minimum over the manifold SO(N) of orthogonal rotations in RN. This is related to the Joint Approximated Diagonalization Eigen-matrices’ (JADE) method [32].The rotation matrixR=(r1,r2,…,rN) composed of the Singular vectors (SVecs) of M(y) (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) ofM(y):η1,N≥η2,N≥⋯≥ηN,N can be expressed as (7)ηi,N=112∑p,q=1NSzipq2+148∑p,q,r=1NKzipqr2;1≤i≤N,which merge all the negentropy terms where the projection of z on the ith SVec: Zi=ri′y 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 Tr[Mk(z)], k≤N (e.g., determinant of M).The vectors corresponding to nonnull SVs span the estimated non-Gaussian subspaceWng-Ed(z) whereas the kernel Ker(M) spanned by the trailing SVecs associated with zero SVs spans the estimated Gaussian subspace Wg-Ed(z) as far as third- and fourth-order cumulants are concerned. Note, however, that the method provides an overestimation of the Gaussian manifold Wg-Ed(z)⊇Wg(z) and an underestimation of the non-Gaussian manifold Wng-Ed(z)⊆Wng(z).To summarize, the first step of source separation is to discard the Gaussian subspace and restrict the analysis to vectors inWng-Ed(z)=Ker[M(z)]⊥ of dimension Nng-Ed. Let us denote the non-Gaussian part of z by zng=Z1,…,ZNng-Ed′.A given partitionzng′=[z1′,z2′,…,zNs′] that leads to separated true non-Gaussian sources, that is, IEd(z(1),z(2),…,z(Ns))=0, makes the invertible matrix M(zng) block-diagonal with Ns blocks. The SVD of M(zng) yields blocks of SVecs: w(i); i=1,…,Ns, spanning the same source’s subspaces (i.e., Span[w(i)]=Span[z(i)]; i=1,…,Ns) (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 inWng-Ed by maximizing the NE sum ∑i=1Nng-EdJEd(Zi) in (6), that is, solving an ICA problem [37, 38] providing independent components (ICs) Zi.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 ofJEd(zng) into positive terms, hereafter called interactivities (ITIs). The decomposition runs through all possible nonnull Nint≡2m-2m-min(m,pmax) subsets u⊆zng with dimension 1≤|u|≤min(m,pmax=4), where (m≡Nng-Ed). Therefore,(8)JEdzng≡∑u⊆zJEd-ITu,JEd-ITu≡112∑i,j,k:ui,uj,uk=uSuijk2+148∑i,j,k,l:ui,uj,uk,ul=uKuijkl2.As an example, if u=U1,U2′ and keeping in mind the hypersymmetry, we get (9)JEd-ITu=1123S1122+3S2212+1484K11122+6K11222+4K22212.The ITIs, for 1D, that is, dim(u)=1 (self-interactivities), coincide with single NEs whereas for dim(u)≥2 they coincide with joint-interactivities, which provide a measure of the non-linear statistical relationships within u contributing to JEd(zng). Using cumulant properties [39], it is straight forward to show that JEd-IT(u)=0 when u has non-null intersections with different independent sources. Therefore, after computing all ITIs, the source identification consists of determining a partition {u(1),u(2),…,u(Ns)} of z composed of subvectors of minimal dimension such that JEd-IT(v)=0 for all v intersecting at least two subvectors u(α) and u(β), α≠β. We note that the NE decomposition is not unique. In fact, the concept of interaction information [40–42] provides a general frame for the decomposition of J(zng) in terms of interaction information, defined for every u⊆zng, which are given by linear combinations of mutual information of subvectors v⊆u. 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 functionf(z) of z-function: E[f(z)] (e.g., the cumulants (2), the HOSVD matrix M(z) (4), the negentropy (5), and its interactivities (8)), through its estimation obtained by sample averages: f(z)¯. 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: Ndof≤Nt. The bias becomes even larger when we deal with high-order cumulants. In fact, the pth order (hypersymmetric) cumulant tensor in N dimensions has CpN+p-1=(N+p-1)!/[(N-1)!p!] independent entries, which scales as a power law of N. This can lead to a high bias of the cumulants of the leading SVs of M(z) 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 N, or high-order statistics p. 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 p.Another strategy that we will also follow here is to evaluate the distribution, moments, and quantiles of the distribution of the estimatorsf(z)¯G, necessary to perform ISA, under the null hypothesis (HG) of multivariate Gaussianity of z. All quantities under HG will hereafter be noted with the subscript “G.” For iid realizations, all the above statistics are scaled as functions of Nt and N, 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 Ndof, taken here as simple red noise AR(1) processes fitted to each z component as g(t+1)=c^1g(t)+(1-c^12)1/2w(t); t=0,…,Nt, where w(t) is a standard Gaussian white noise process and c^1 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 hypothesisHG. The finite-sample estimator of JEd, hereafter denoted as J~Ed (all quantities based on sample statistics are overlined with a tilde) and others based on truncated NE approximations are conservative tests of HG (or less powerful tests of its negation: ~HG), because they reject non-Gaussianity (Type II error) resulting from discarded cumulants of order p>pmax=4.Now, we discuss the behavior of the estimatorJ~Ed:G of JEd:G (under HG) computed from Nt-sized samples of a vector of N independent scalars. The pdf of J~Ed:G is positively skewed, depending on the joint distribution of S~G and K~G, which are in general not independent since they satisfy a number of inequality constraints (e.g., K1111-S1112+2≥0). The bias of J~Ed:G is given by(10)EJ~Ed:G=1NdofPJEdN+O1Ndof2,1NdofPJEdN≡1121NdofN3+3N2+2N+1481NdofN4+6N3+11N2+6N,where the first and second square-bracket terms of (10) represent, respectively, scaled versions of E[S~GF2] and E[K~GF2]. The bias (10) tends asymptotically to PJEd(N)/Ndof. It also provides a lower-bound of the positive J~Ed 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 3N(N-1) cumulant biases of type E[(YiYiYj)¯2]=2Ndof-1+O(Ndof-2); i≠j in E[S¯GF2] based on the moments of the standard Gaussian: E(U2k)=(2k)!/(2kk!). Figure 1 shows the dependence on Nt and c^1 for N=4, of the ratio(11)fREdNt,c^1,N≡EJ~Ed:GNtPJEdNbetween the actual bias and the asymptotic bias of J~Ed:G (at c^1=0), obtained from 1000 Monte-Carlo replicates. Moreover, the dependence of fREd on N is quite negligible with relative deviations less than 5% (not shown).Figure 1
(a) Contour plot offREd as a function of Nt, c^1 for a fixed dimension N=4. (b) Ratio fREd-Psig for the confidence level psig=0.95, as a function of the dimension N and the lag-one autocorrelation c^1 for a sample size Nt=1000. The ratio increases with decreasing N or increasing c^1.
(a)
(b)From the analysis of Figure1(a) and as expected, the ratio fREd remains close to 1 when c^1 is small and Nt is large. For Nt≥500, the ratio has reached the asymptotic value, where it becomes dependent only on c^1 with Ndof=Nt/fREd(∞,c^1,N)<Nt, where fREd is interpreted as the “average interval between independent realizations.” For small sample sizes, fREd<fREd(∞,c^1,N), due to a negative higher-order term in (10).The variance ofJ~Ed:G involves covariances between generic cumulants, for example, cm(u1) and cm(u2), of subvectors u1,u2⊆z. Those covariances satisfy cov[cm(u1),cm(u2)]=kNt-d; k≥0 with d=2 if u1=u2 and d=3 if u1⊃u1∩u2≠⌀ and d=∞ otherwise.As for the bias, the quantilespsig∈[0,1] of J~Ed:G hereby noted Qpsig(J~Ed:G) scale asymptotically as Nt-1. Figure 1(b) shows the dependence on N and c^1 for Nt=1000, of the ratio fREd-Psig≡Qpsig(J~Ed:G)[Nt/PJEd(N)]. This figure provides a kind of a rule of thumb thresholds of J~Ed for rejecting HG in the asymptotic regime.


### 2.2.2. Estimation of the Gaussian Subspace by HOSVD

Following a rejection ofHG, the next step is to estimate the Gaussian subspace Wg and its dimension Ng. With finite samples, one computes the sample-based matrix M~(z) of M(z) (4). However, Wg can no longer be estimated by Ker(M~) (see Section 2.1) but through the subspace spanned by SVecs associated with the trailing (not significant) SVs.UnderHG, the (estimated) SVs η~i,N:G, (see (7)), with quantiles Qpsig(η~i,N:G), i=1,…,N, are asymptotically O(Ndof-1) and the SVs become undistinguishable for large Ndof [36]. However, the SVs are not independent and consequently their marginal quantiles cannot be estimated independently for each i. To overcome this difficulty, a sequence of statistical tests is applied, as discussed in the next subsection.In summary, we first consider the vectorsz=Z1,…,ZN′=R′y and zs≡Zs,1,…,Zs,N′=Rs′ys, where ys is a generic surrogate of y and R and Rs are matrices filled by the ranked SVecs of M~(y) and M~(ys), respectively. For each i in {1,…,N}, we construct the vector bi≡(Z1,…,Zi-1,Zs,i,…,Zs,N), where the trailing i,…,N components are swapped by the associated ranked surrogates. Then we compute their NE contributions M~ii(bi) (5) and finally the quantile Qpsig[M~ii(bi)] evaluated from a large ensemble of surrogates. The test then proceeds backward, starting from the smallest SV (i=N). Therefore, if M~kk(z)<Qpsig[M~kk(bk)], for all k≥i, then the null hypothesis HG,i (Zkundistinguishable from Gaussian surrogates Zs,k for all k≥i) cannot be rejected. The smallest i for which HG,i is not rejected determines the estimated Gaussian subspace Wg-Ed with dimension Ng-Ed (spanned by the trailing Ng-Ed SVecs of R), and Wng-Ed≡Wg-Ed⊥ (spanned by the leading Nng-Ed≡N-Ng-Ed SVecs of R).


### 2.2.3. Estimation of Multivariate Sources

The next step is to identify the non-Gaussian sources within the vector of non-Gaussian componentszng≡Z1,…,ZNng-Ed′ spanning Wng-Ed. Since zng is a subvector of z, we easily get J~Ed(zng)≤J~Ed(z).We test the separability ofWng-Ed into scalar and/or subvector sources by looking for a new orthogonal rotation zISA≡RISA′zng=Z1,ISA,…,ZNng-Ed,ISA′=z1,ISA′,z2,ISA′,…,zNs,ISA′′ composed of a maximal number Ns≥2 of subvector candidate sources z(i),ISA, i=1,…,Ns, such that the joint pdf satisfies(12)ρzISA=1+εISA∏k=1Nsρzk,ISA,where the norm |εISA| of the pdf separation error εISA is below the threshold of rejection of the null hypothesis of nonseparability (Ns=1). There is an a priori source configuration indeterminacy. For instance, for Nng-Ed=4, 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 is(13)FISAR≡∑i=1NsJ~Edzi,ISA=J~Edzng-I~Edz1,ISA,z2,ISA,…,zNs,ISA,where I~Ed 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 (14)RICA≡argmaxRISA∑i=1NgJ~EdZi,ISA,providing the vector of best scalar sources: zICA=RICA′zng=Z1,ICA,…,ZNNg-Ed,ICA′. The algorithm to get RICA is explained in Appendix A. Finally, we build multivariate sources forming appropriate groups of zICA components, according to the conjecture of [13], and as below.We start with the NE decomposition intoNint interactivities (8) as(15)J~Edzng=J~EdzICA=∑u⊆zICAJ~Ed-ITu.To compare and account for the statistical significance of u-interactivities in (15) with different dimensions dim⁡(u), we introduce scaled interactivities J~Ed-ITsig(u) through its comparison with the quantiles of the NE J~Ed-IT(uG) of the ranked (by HOSVD) Gaussian surrogates uG of u. It is is defined as(16)J~Ed-ITsigu=J~Ed-ITwQpsigJ~Ed-ITuG;w=uifdimu=1uGsiifdimu>1.When dim(u)>1 in (16), w is the rotation of uGs towards isotropy (e.g., by Gram Schmidt orthogonalization) with uGs being composed of the standard Gaussianized components of u. Since uG has Gaussian marginals, this procedure leads to a fairer and unbiased test when u has non-Gaussian pdf marginals. This is also supported by the fact that the mutual information of u is I(u)=I(uGs)=I(uGsi)+[-1/2log⁡(det⁡CGs)]+∑k=1NJ(uGsi,k)≈I(uGsi) where the second term is the Gaussian MI component [47, 48], which depends on the correlation matrix CGs of uGs and the third term is the sum of the NEs of the components of uGsi. The above measure is therefore resistant to single outliers of the marginal pdfs.The differentNint scaled interactivities (16) are then sorted by decreasing order, retaining uniquely the leading ones verifying J~Ed-ITsig(u)≥Jsig where Jsig 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, ∑uJ~Ed-ITsig(u)<J~Ed(zng)-E[J~Ed(zng)], which will be adopted in practice.


## 2.2.1. Statistical Tests of Multivariate Non-Gaussianity

The above source separation method depends crucially on the approximation of the expectation of any functionf(z) of z-function: E[f(z)] (e.g., the cumulants (2), the HOSVD matrix M(z) (4), the negentropy (5), and its interactivities (8)), through its estimation obtained by sample averages: f(z)¯. 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: Ndof≤Nt. The bias becomes even larger when we deal with high-order cumulants. In fact, the pth order (hypersymmetric) cumulant tensor in N dimensions has CpN+p-1=(N+p-1)!/[(N-1)!p!] independent entries, which scales as a power law of N. This can lead to a high bias of the cumulants of the leading SVs of M(z) 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 N, or high-order statistics p. 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 p.Another strategy that we will also follow here is to evaluate the distribution, moments, and quantiles of the distribution of the estimatorsf(z)¯G, necessary to perform ISA, under the null hypothesis (HG) of multivariate Gaussianity of z. All quantities under HG will hereafter be noted with the subscript “G.” For iid realizations, all the above statistics are scaled as functions of Nt and N, 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 Ndof, taken here as simple red noise AR(1) processes fitted to each z component as g(t+1)=c^1g(t)+(1-c^12)1/2w(t); t=0,…,Nt, where w(t) is a standard Gaussian white noise process and c^1 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 hypothesisHG. The finite-sample estimator of JEd, hereafter denoted as J~Ed (all quantities based on sample statistics are overlined with a tilde) and others based on truncated NE approximations are conservative tests of HG (or less powerful tests of its negation: ~HG), because they reject non-Gaussianity (Type II error) resulting from discarded cumulants of order p>pmax=4.Now, we discuss the behavior of the estimatorJ~Ed:G of JEd:G (under HG) computed from Nt-sized samples of a vector of N independent scalars. The pdf of J~Ed:G is positively skewed, depending on the joint distribution of S~G and K~G, which are in general not independent since they satisfy a number of inequality constraints (e.g., K1111-S1112+2≥0). The bias of J~Ed:G is given by(10)EJ~Ed:G=1NdofPJEdN+O1Ndof2,1NdofPJEdN≡1121NdofN3+3N2+2N+1481NdofN4+6N3+11N2+6N,where the first and second square-bracket terms of (10) represent, respectively, scaled versions of E[S~GF2] and E[K~GF2]. The bias (10) tends asymptotically to PJEd(N)/Ndof. It also provides a lower-bound of the positive J~Ed 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 3N(N-1) cumulant biases of type E[(YiYiYj)¯2]=2Ndof-1+O(Ndof-2); i≠j in E[S¯GF2] based on the moments of the standard Gaussian: E(U2k)=(2k)!/(2kk!). Figure 1 shows the dependence on Nt and c^1 for N=4, of the ratio(11)fREdNt,c^1,N≡EJ~Ed:GNtPJEdNbetween the actual bias and the asymptotic bias of J~Ed:G (at c^1=0), obtained from 1000 Monte-Carlo replicates. Moreover, the dependence of fREd on N is quite negligible with relative deviations less than 5% (not shown).Figure 1
(a) Contour plot offREd as a function of Nt, c^1 for a fixed dimension N=4. (b) Ratio fREd-Psig for the confidence level psig=0.95, as a function of the dimension N and the lag-one autocorrelation c^1 for a sample size Nt=1000. The ratio increases with decreasing N or increasing c^1.
(a)
(b)From the analysis of Figure1(a) and as expected, the ratio fREd remains close to 1 when c^1 is small and Nt is large. For Nt≥500, the ratio has reached the asymptotic value, where it becomes dependent only on c^1 with Ndof=Nt/fREd(∞,c^1,N)<Nt, where fREd is interpreted as the “average interval between independent realizations.” For small sample sizes, fREd<fREd(∞,c^1,N), due to a negative higher-order term in (10).The variance ofJ~Ed:G involves covariances between generic cumulants, for example, cm(u1) and cm(u2), of subvectors u1,u2⊆z. Those covariances satisfy cov[cm(u1),cm(u2)]=kNt-d; k≥0 with d=2 if u1=u2 and d=3 if u1⊃u1∩u2≠⌀ and d=∞ otherwise.As for the bias, the quantilespsig∈[0,1] of J~Ed:G hereby noted Qpsig(J~Ed:G) scale asymptotically as Nt-1. Figure 1(b) shows the dependence on N and c^1 for Nt=1000, of the ratio fREd-Psig≡Qpsig(J~Ed:G)[Nt/PJEd(N)]. This figure provides a kind of a rule of thumb thresholds of J~Ed for rejecting HG in the asymptotic regime.


## 2.2.2. Estimation of the Gaussian Subspace by HOSVD

Following a rejection ofHG, the next step is to estimate the Gaussian subspace Wg and its dimension Ng. With finite samples, one computes the sample-based matrix M~(z) of M(z) (4). However, Wg can no longer be estimated by Ker(M~) (see Section 2.1) but through the subspace spanned by SVecs associated with the trailing (not significant) SVs.UnderHG, the (estimated) SVs η~i,N:G, (see (7)), with quantiles Qpsig(η~i,N:G), i=1,…,N, are asymptotically O(Ndof-1) and the SVs become undistinguishable for large Ndof [36]. However, the SVs are not independent and consequently their marginal quantiles cannot be estimated independently for each i. To overcome this difficulty, a sequence of statistical tests is applied, as discussed in the next subsection.In summary, we first consider the vectorsz=Z1,…,ZN′=R′y and zs≡Zs,1,…,Zs,N′=Rs′ys, where ys is a generic surrogate of y and R and Rs are matrices filled by the ranked SVecs of M~(y) and M~(ys), respectively. For each i in {1,…,N}, we construct the vector bi≡(Z1,…,Zi-1,Zs,i,…,Zs,N), where the trailing i,…,N components are swapped by the associated ranked surrogates. Then we compute their NE contributions M~ii(bi) (5) and finally the quantile Qpsig[M~ii(bi)] evaluated from a large ensemble of surrogates. The test then proceeds backward, starting from the smallest SV (i=N). Therefore, if M~kk(z)<Qpsig[M~kk(bk)], for all k≥i, then the null hypothesis HG,i (Zkundistinguishable from Gaussian surrogates Zs,k for all k≥i) cannot be rejected. The smallest i for which HG,i is not rejected determines the estimated Gaussian subspace Wg-Ed with dimension Ng-Ed (spanned by the trailing Ng-Ed SVecs of R), and Wng-Ed≡Wg-Ed⊥ (spanned by the leading Nng-Ed≡N-Ng-Ed SVecs of R).


## 2.2.3. Estimation of Multivariate Sources

The next step is to identify the non-Gaussian sources within the vector of non-Gaussian componentszng≡Z1,…,ZNng-Ed′ spanning Wng-Ed. Since zng is a subvector of z, we easily get J~Ed(zng)≤J~Ed(z).We test the separability ofWng-Ed into scalar and/or subvector sources by looking for a new orthogonal rotation zISA≡RISA′zng=Z1,ISA,…,ZNng-Ed,ISA′=z1,ISA′,z2,ISA′,…,zNs,ISA′′ composed of a maximal number Ns≥2 of subvector candidate sources z(i),ISA, i=1,…,Ns, such that the joint pdf satisfies(12)ρzISA=1+εISA∏k=1Nsρzk,ISA,where the norm |εISA| of the pdf separation error εISA is below the threshold of rejection of the null hypothesis of nonseparability (Ns=1). There is an a priori source configuration indeterminacy. For instance, for Nng-Ed=4, 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 is(13)FISAR≡∑i=1NsJ~Edzi,ISA=J~Edzng-I~Edz1,ISA,z2,ISA,…,zNs,ISA,where I~Ed 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 (14)RICA≡argmaxRISA∑i=1NgJ~EdZi,ISA,providing the vector of best scalar sources: zICA=RICA′zng=Z1,ICA,…,ZNNg-Ed,ICA′. The algorithm to get RICA is explained in Appendix A. Finally, we build multivariate sources forming appropriate groups of zICA components, according to the conjecture of [13], and as below.We start with the NE decomposition intoNint interactivities (8) as(15)J~Edzng=J~EdzICA=∑u⊆zICAJ~Ed-ITu.To compare and account for the statistical significance of u-interactivities in (15) with different dimensions dim⁡(u), we introduce scaled interactivities J~Ed-ITsig(u) through its comparison with the quantiles of the NE J~Ed-IT(uG) of the ranked (by HOSVD) Gaussian surrogates uG of u. It is is defined as(16)J~Ed-ITsigu=J~Ed-ITwQpsigJ~Ed-ITuG;w=uifdimu=1uGsiifdimu>1.When dim(u)>1 in (16), w is the rotation of uGs towards isotropy (e.g., by Gram Schmidt orthogonalization) with uGs being composed of the standard Gaussianized components of u. Since uG has Gaussian marginals, this procedure leads to a fairer and unbiased test when u has non-Gaussian pdf marginals. This is also supported by the fact that the mutual information of u is I(u)=I(uGs)=I(uGsi)+[-1/2log⁡(det⁡CGs)]+∑k=1NJ(uGsi,k)≈I(uGsi) where the second term is the Gaussian MI component [47, 48], which depends on the correlation matrix CGs of uGs and the third term is the sum of the NEs of the components of uGsi. The above measure is therefore resistant to single outliers of the marginal pdfs.The differentNint scaled interactivities (16) are then sorted by decreasing order, retaining uniquely the leading ones verifying J~Ed-ITsig(u)≥Jsig where Jsig 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, ∑uJ~Ed-ITsig(u)<J~Ed(zng)-E[J~Ed(zng)], 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 (https://www.ncdc.noaa.gov/oa/climate/research/sst/ersstv3.php), from the International Comprehensive Ocean–Atmosphere Data Set (ICOADS) release 2.4. The SST data are on 2°×2° 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 Nt=1224 monthly SST anomalies (SSTAs), on Np=8714 grid points. A EOF analysis is then applied [51] and we only use the space spanned by the leading Nmax=12 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 σTotal2=1709.4 (°C)2, the local variance at point P is σP2, with an average σLocal2≡σTotal2/[NPclat]~(0.50)2 (°C)2 (clat = average of cos⁡(latitude)) and we let the fractional local variance σ^P2≡σP2/σLocal2. 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 c^1 based on an AR(1) model.Figure 2
(a) % of explained variance confidence interval and (b) lag-one autocorrelation of the AR(1) fitting model of the leading 12 PCs.
(a)
(b)Figure3 shows the leading 12 EOFs which are described, with references in Table 1. Based on c^1, an estimate of the “average interval between independent realizations” is Δind=Nt/Ndof=fREd(Nt,c^1,N), (Figure 1(a)) which is not too different from (-1.)/ln⁡(c^1) (e.g., [53]).Table 1
EOF description with the main SSTA modes, along with corresponding available indices from literature and references.
EOF rank
Description with the main SSTA modes
Available index
PC/Index cor.
Refs.
Notes
1
El Niño (ENSO)
MEI
0.98
[56]
2
Atlantic Multi-Decadal Oscillation (AMO)
AMO
0.80
[57]
3
Pacific Decadal Oscillation (PDO): negative phase
PDO
-
0.56
[58, 59]
4
Northern Pacific Gyre Oscillation (NPGO)
NPGO
0.53
[60]
El Niño Modoki (EM): negative phase
EMI
-
0.42
[61]
Also known as Central Pacific or noncanonical El Niño
5
El Niño Modoki (EM): positive phase
EMI
0.57
[61]
Subtropical Indian Ocean Dipole: negative phase
SIOD
[62, 63]
First EOF mode of the SSTA in the Indian Ocean ranging in the 15°S–45°S band
6
Atlantic Niño: negative phase
ATL3
-
0.49
[64, 65]
South Atlantic Subtropical Dipole
SASD
[66]
First EOF coupled mode of the SSTA/wind in the South Atlantic
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
Zonally oriented dipole in South Pacific
[67, 68]
Second EOF of the South Pacific SSTA
7
Meridional sequence of dipoles crossing the North-South Pacific
South Pacific Ocean Dipole (SPOD)
SPOD
[67, 69]
First EOF of the South Pacific SSTA
8
North Atlantic tripole (NAT)
NAT
[70]
Tropical Atlantic SST dipole
[68]
9
South Pacific Ocean Dipole (SPOD)
SPOD
[67, 69]
First EOF of the South Pacific SSTA
10
Zonally oriented dipole in South Atlantic
[68]
Zonally oriented dipole in South Pacific
[67, 68]
Second EOF of the South Pacific SSTA
11
Meridional 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
12
Tripolar pattern in North Pacific
North Atlantic tripole (NAT)
NAT
[70]Table 2
% of explained variance, skewness, and kurtosis of the 12 leading PCs and the squared projections ofWng-Ed onto the individual whitened PCs for the analysis over N=11 PCs (see the end of Section 3.2).
PC1
PC2
PC3
PC4
PC5
PC6
PC7
PC8
PC9
PC10
PC11
PC12
%Var
17.6
5.3
4.6
3.8
3.3
2.7
2.5
2.2
2.2
2.0
1.9
1.8
Skew.
0.24
0.03
0.01
0.28
-
0.16
0.25
-
0.08
0.23
0.14
0.16
0.29
0.13
Kurt.
-
0.01
-
0.64
-
0.33
0.19
0.02
-
0.02
-
0.06
0.17
-
0.04
0.55
0.14
-
0.16
SqPro
0.65
0.53
0.16
0.68
0.50
0.56
0.29
0.05
0.37
0.81
0.39Figure 3
Maps of the 12 variance-leading, spatially normalized EOFs of the monthly SSTAs field, computed in the period 1910–2011. Ranked EOFs are read from left to right and top to bottom.Figure2(b) shows that Δind 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 c^1 values, that is, 0.86 (see Section 3.3).


### 3.2. Negentropy Distribution and Non-Gaussian Subspace Estimation

Non-Gaussianity and negentropy of a fieldx^ 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 y, through linear operators. Let us consider, for example, the local-basis coskewness tensor S(x^). 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 (P,Q,R):(17)Sx^F2≡∑P,Q,R=1NPSx^PQR2σ^P2σ^Q2σ^R2a^Pa^Qa^R=∑i,j,k=1NrSyijk2λ^iλ^jλ^k≃SyF2Nr3∑i,j,k=1Nrλ^iλ^jλ^k,where a^P stand for the relative area at grid point P and σ^P2 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 NEJ~Ed(y) can be decomposed using M~(y) (5). For each N=1,…,Nmax=12, we denote by y(N)=(Y1,…,YN) the set of N leading whitened PCs. Figure 4 shows the diagonal terms (M~kk[y(N)]), Figure 4(a) and the total NE J~Ed[y(N)]=Tr(M~[y(N)]), Figure 4(b). Note that, for a given N, the smaller values of M~kk[y(N)] are either the least non-Gaussian PCs or the less interactive PCs with the remaining PCs (i≠k), providing therefore good candidates for the Gaussian subspace of y(N). For example, with N=12, we have k=2,7,8,9, 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 M~22[y(N)] of the PC2 due to its small degrees of freedom.Figure 4
(a) NE components parcelM~kk[y(N)] (k in abscissas) for each N=1,…,12 (marked in grey shaded). (b) Total NE J~Ed[y(N)] for N=1,…,12. Values rejecting the null hypothesis of multi-Gaussianity at confidence levels <70%, 70%–90%, 90%–95%, and >95% are marked, respectively, by “+,” open circle, small solid circle, and large solid circle.
(a)
(b)BiasesE{J~Ed[y(N)]}-JEd[y(N)] may be approximated as E(J~Ed:G)=fREd(Nt,c^1,N)PJEd(N)/Nt (see (11)) with c^1≈0.86, leading to fREd(Nt=1224,c^1=0.86)~3.3 (see Figure 1(a)). For instance, for N=12, J~Ed=3.1, 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.Figure4(a) shows that, for all N, the largest NE component is M~11[y(N)], 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:M~kk[y(N)]-M~kk[y(N-1)] (Figure 4(a)), one may detect which interactivities mostly contribute to J~Ed[y(N)]. The largest differences are for (k,N) pairs: {(1,3),(1,5),(1,10),(1,11), (3,11),(4,10),(4,11),(5,10),(5,11),(6,11)} in consistency with the high confidence level of M~kk[y(11)] for k=1,3–6,10,11. Those jumps are often due to particular joint cumulants in the 2D subspaces spanned by (PCk,PCN), responsible for particular types of joint non-Gaussianity, for example, scatter-plots resulting for instance from the occurrence of joint extremes. For instance, Yk2Yi¯>0, is favored for the sign combinations (-,+) and (+,+) of (Yk,Yi). As an illustration, Figures 5(a) and 5(b) show, respectively, scatter-plots of (Y1g,Y3g) and (Y1g,Y11g) (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 3 of [50]) during which a positive PC11 and a negative PDO phase have occurred.Figure 5
Scatter-plots of PC1 versus PC3 (a) and PC1 versus PC11 (b). The marginal pdfs have been Gaussianized.
(a)
(b)SinceM~12,12[y(12)] 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 Hg test to limit N [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 vectory(11) and the associated vector z=Z1,…,Z11′, yielding J~Ed[z]=2.55, and a bias of 1.73, leaving ~0.82 for “true non-Gaussianity.” Figure 6 shows the NE components (or SVs) M~kk(z) (k=1,…,N), along with the threshold quantiles 90% and 95% of Qpsig[M~kk(bk)] obtained from an ensemble of 1000 Gaussian surrogates. For comparison, we also show the SVs M~kk(zs) (open squares) of surrogate realizations zs of z. The SVs of zs all lie below the above thresholds, supporting the consistency of the Gaussianity test. Figure 6 also suggests that the estimated Gaussian subspace Wg-Ed of the SSTA is spanned by the last 6 SVecs (with a confidence level psig=95%) and so Ng-Ed=6 and, Nng-Ed=5 for the non-Gaussian subspace Wng-Ed. The total NE restricted to zng≡Z1,…,Z5′ is J~Ed[zng]=1.11 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 M~kk(zng) (open black circles), with their low values compared to M~kk(z), k=1,…,5, are also shown in Figure 6.Figure 6
Ranked SVsM~kk(z), k=1,…,N=11, contributing to NE, both from real SSTA data (solid squares) and from Gaussian surrogate realizations (open squares) and 90% (red), 95% (blue) quantiles of Qpsig[M~kk(bk)] obtained from Gaussian surrogates after rank k, to be compared with M~kk(z). The trailing 6 SVs is not significant at psig=95%, leading to Nng-Ed=5. The NE components restricted to the non-Gaussian subspace are marked as open circles.The squared projections ofWng-Ed onto each of the 11 whitened PCs are given by the inner product (Appendix B) Wng-Ed·ei=∑k=15Ri,k2, for i=1,…,N=11, (Table 2), where ei is the (unitary) ith vector. By definition, those values add up to Nng-Ed=5. 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 ∑i=111λ^i∑k=15Ri,k2 = 24.7% and 23.3%, for the non-Gaussian and Gaussian subspaces.


### 3.3. Application to Synthetic Non-Gaussian Data

We generateNmax=12 AR(1) independent series Gk(t), t=1,…,Nt, k=1,…,Nmax. The whitened and uncorrelated components of the non-Gaussian vector y′=(Y1,…,YN) are computed following(18)Yk=Gk;k=1,…,N≤Nmax,k≠2,5Y2=σ2-1/2G12-1+1-σ21/2G2;Y5=σG3G4+1-σ21/2G5with σ2∈[0,1] measuring the “strength” of non-Gaussianity. So y is composed of Ng=N-5 Gaussian scalar sources spanning Wg and Nng=5 non-Gaussian variables spanning Wng. The latter is composed of one dyad y′(1)=(Y1,Y2) and a non-Gaussian triad y′(2)=(Y3,Y4,Y5).The next step consists of testing the null hypothesisHG of global Gaussianity (i.e., σ2=0) and evaluating J~Ed(y) by varying σ2 in [0,1] and Ng in {0,…,7}.Figure7(a) shows the Monte-Carlo average E[J~Ed(y)] obtained from 1000 surrogates. As expected, the mean NE increases with σ2 (nonlinearity effect) and Ng (dimensionality effect), though at a much smaller rate. The bias E[J~Ed:G(y)] is the value reached at σ2=0 in agreement with (11).Figure 7
(a) Monte-Carlo average of the sample negentropy estimatorJ~Ed(y) as a function of σ2 and Ng. (b) Contour plot of the fraction F of times for which the null hypothesis HG is rejected at psig=99% (1% significance level) as a function of σ2and Ng. The contour F=0.99 can mark σmin2 at the 99% confidence level.
(a)
(b)


#### 3.3.1. Gaussian Subspace Recovery

We investigate the strength ofJ~Ed(y) as a function of σ2 and Ng. The null hypothesis HG is rejected at the confidence level psig if J~Ed(y)>[PJEd(N)Nt-1]fREd-Psig using fREd-Psig (Figure 1(b)). The Monte-Carlo fraction F of HG rejection is shown in Figure 7(b) for psig=99%. The contour line at F=0.99 provides the graph σmin2(Ng) of σ2 for the recovery of the non-Gaussian sources. As expected, σmin2 increases with Ng, starting at σmin2=20% and at Ng=0 and reaching σmin2~45% at Ng=6. This means that the larger N, the stronger (through σ2) the non-Gaussian embedded sources. Moreover, σmin2 diminishes under a less conservative test (lower psig) 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., psig=90%) [47].


#### 3.3.2. ICA and ISA

Here, we apply ICA and ISA forN=11 (Ng=6) and compute the skill of source separation in two cases: (a) σ2=20%, high noise, and (b) σ2=40%, intermediate noise, yielding, respectively, total NE of 1.79 and 2.57.Figures8(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: Nng-Ed=5 with psig=90% 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 Wng-Ed·Wng (compared to 5 in the ideal recovery). Those products are, respectively, 3.99 and 4.49 for cases (a) and (b). Lower values of σ2 lead to a subestimation of Nng and to a worse estimation of Wng.Figure 8
Same as Figure6 but for the for synthetic data with σ2=20% (a) and σ2=40% (b).
(a)
(b)Now, by applying ICA (see Section2.2.3 and Appendix A) we get ICs making the vector zICA=(RICA′R′)y=Z1,ICA,…,ZNNg-Ed,ICA′. Tables 3 and 4 give, respectively, for cases (a) and (b), a summary of ICA with the sorted maximized self-NEs J~Ed(Zi,ICA), the NE components M~ii(zICA), and the IC vector loadings in RRICA. 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 zEd-dyad is spanned by IC1 and IC4 whereas the estimated triad zEd-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 σ2. 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).Table 3
Contributions to the ICA contrast function (2nd column). NE components (3rd column) and IC vector loadings on the basis of original variables (18) for the high noise case. Components in which ICs are most projected are in bold font.
IC rank
J
~
E
d
(
Z
i
,
I
C
A
)
M
~
i
i
(
z
I
C
A
)
Y
1 (dyad)
Y
2 (dyad)
Y
3 (triad)
Y
4 (triad)
Y
5 (triad)
IC1:i=1
0.127
0.162
0.76
-
0.58
−0.02
0.01
0.16
IC2:i=2
0.066
0.125
−0.07
0.07
0.47
-
0.56
0.56
IC3:i=3
0.042
0.102
−0.01
−0.08
-
0.38
-
0.68
-
0.42
IC4:i=4
0.014
0.053
-
0.41
-
0.52
0.21
0.29
0.01
IC5:i=5
0.011
0.063
−0.03
0.02
-
0.68
0.07
0.57
Total
0.260
0.506Table 4
Table3 but for the case of intermediate noise (b).
IC rank
J
~
E
d
(
Z
i
,
I
C
A
)
M
~
i
i
(
z
I
C
A
)
Y
1 (dyad)
Y
2 (dyad)
Y
3 (triad)
Y
4 (triad)
Y
5 (triad)
IC1:i=1
0.410
0.454
0.70
-
0.70
−0.03
0.00
0.08
IC2:i=2
0.256
0.372
0.00
0.05
0.45
0.65
0.58
IC3:i=3
0.109
0.257
0.03
−0.03
-
0.30
0.72
-
0.55
IC4:i=4
0.033
0.078
0.59
0.59
−0.17
−0.07
0.00
IC5:i=5
0.031
0.114
−0.04
0.04
-
0.76
0.11
0.54
Total
0.840
1.277The final step consists of merging the ICs as outlined in Section2.2.3. Figure 9 presents the results for both cases, where the self-interactivities (16) are the most relevant ITIs due to their maximization by ICA (except IC5 in case (a)). The cumulated NE reaches the unbiased NE value (dashed line) near the threshold of the ITI acceptance (scaled-ITI > 1). In case (a) (Figure 9(a)), the accepted ITIs are for ICs 1–5, and dyadic links (IC3, IC2), (IC5, IC3) yielding the emergence of the true triad (IC2, IC3, and IC5) (see Table 3). The marginally accepted ITIs holds for (IC4, IC3, and IC2) (false triad) and (IC1, IC4) (true dyad). Excluded ITIs remain at the level of the NE bias. For smaller values of σ2 (not shown), the recovered sources may be mixed, partly restored, or even not detected.Figure 9
Sorted scaled-interactivities by decreasing order (solid circles), ITIs (solid squares) with cumulated ITIs (open crossed circles) and unbiased total NE (horizontal dashed line). The sorted IC groups are labeled (e.g., 32 means group of IC2 and IC3). The 95% confidence levels of the selected ITIs are those for which the cumulated ITIs are below the dashed line. (a) Case (a): high noise and (b) case (b): intermediate noise.
(a)
(b)For case (b) (Figure9(b)) the source recovery improves. In fact, the significant accepted ITIs are for ICs 1–5, the links for (IC5, IC3), (IC5, IC2), (IC3, IC2) with the true emerging triad (IC2, IC3, and IC5), and for the true dyad (IC4, IC1). The sum of ITIs associated with ICA and ISA is, respectively, 0.84 and 1.13, showing clearly the advantage of ISA in this case.


### 3.4. Application of ICA and ISA to SSTA Data

Using the estimated non-Gaussian subspace (Section3.2), we apply here ICA to the leading NICA=5 SVecs of the SSTA dataset (see Figure 6). Table 5 gives a summary of the application with the characteristics of the ranked ICs Zi,ICA, i=1,…,5, by decreasing self-NEs. They explain, respectively, 7.1%, 5.5%, 4.3%, 4.1%, and 3.7% of the total variance.Table 5
Sorted ICsi=1,…,5, for SSTA data showing self-NEs contributing to the ICA contrast function (2nd column), NE components (3rd column), and loading factors (RRICA)k,i (i=1,…,5), of the leading ICs on the PC basis. Loading factors, sorted by absolute value, explaining up to 85% of the unitary norm are in bold font. The last two columns give the skewness and kurtosis of ICs.
IC rank
J
~
Ed
(
Z
i
,
I
C
A
)
M
~
i
i
(
z
I
C
A
)
PC1
PC2
PC3
PC4
PC5
PC6
PC7
PC8
PC9
PC10
PC11
Sk.
Kurt.
IC1:i=1
0.694
0.747
0.53
0.05
0.34
0.40
−0.17
0.24
0.06
0.12
−0.10
0.18
0.52
1.48
4.95
IC2:i=2
0.099
0.137
0.44
0.13
−0.14
0.02
0.03
-
0.63
-
0.35
−0.10
-
0.44
0.13
−0.15
−0.62
1.79
IC3:i=3
0.042
0.091
0.24
-
0.44
0.11
0.35
0.55
-
0.27
0.24
0.03
0.28
-
0.29
−0.12
−0.44
1.12
IC4:i=4
0.042
0.074
−0.26
0.35
−0.07
0.47
0.36
0.12
0.28
0.10
-
0.29
0.42
-
0.27
−0.24
1.32
IC5:i=5
0.016
0.064
−0.21
0.44
0.06
0.42
−0.18
−0.12
−0.14
0.10
−0.04
-
0.70
0.03
−0.32
0.58
Total
0.893
1.113Most of the NE components come from the self-NEs (an effect of ICA), through their high skewness and positive kurtosis or leptokurtosis (Table5), leading to fat long pdf tails. Extreme values of about 3–5 standard deviations are often reached as seen by the IC time-series (Figure 10). Under those conditions, ICs are basically aggregative and cooperative combinations of single and joint PC extremes. Quasi-independence and noncorrelation between ICs means that occurrence of extremes among different ICs tends to be asynchronous (see Figure 10) and therefore extreme phenomena are independent.Figure 10
IC time-series for the period 1910–2011.Figure11 shows the correlation maps between the ICs and SSTAs. An alternative way to interpret ICs is through the patterns displayed in the corresponding normalized loading maps (Figure 12) obtained from combinations of EOFs (Figure 3). Correlation and loading map values at the grid point P for the ith IC are weighted by the PCs standard deviation and its inverse, respectively, ∑k=1N=11Dk,Pek(P)(RRICA)k,i where Dk,P is λ^k1/2/σ^P (Figure 11) and λ^k-1/2 (Figure 12), respectively, and ek(P) is the kth EOF at P (Figure 3). By definition, ICs are spatial inner products between correspondent normalized loading maps (Figure 12) and the SSTA field. Consequently, maximal positive (negative) ICs on a given month are favored by highly positive (negative) pattern correlations between their loading maps and the SSTA field. Note, however, that the patterns shown in Figures 11 and 12 are not orthogonal. On loading maps, the weighting gives rise to the impact on ICs of PCs of small explained variance. Because of this weighting, the correlation maps (Figure 11) often display much broader and large-scale patterns than loading maps (Figure 12) [15], though both share the same main features when the range of PC variances contributing to a given IC is not too wide, which is the case here.Figure 11
Correlation maps between ICs and SSTA field.Figure 12
Spatially normalized loading maps of the 5 sorted ICs.The leading IC1 explains about 2/3 of the total NE with the remaining 4 being almost evenly distributed (Table5). Further experiments show that IC1 is quite robust with the use of NICA<5 leading SVecs with the remaining ICs being slight mixtures of those obtained with NICA=5. The above mixture of PCs suggests the existence of a global nonseparable source thanks to the nonlinear interaction across scales. This strongly suggests that ICs, obtained as blind source separation, are convenient tools for identifying nonlinear variability. Dynamical constraints may be included explicitly, for example, via nonlinear dynamical modes [71, 72].The leading IC1 reaches values of +5 standard deviations in the years 1983 and 1997 of strong El Niños, during which PC1, PC3, PC4, and PC11 are synchronised (see also Figures5(a) and 5(b) and Table 5). The coherence between those four PCs may sometimes be partial.The correlation pattern for IC1 (Figure11) has an El Niño pattern shape due to the high PC1 amplitude (Table 5) and PC1 variance whereas the loading map (Figure 12) is much correlated with EOF11, showing more local features. Therefore, from Figure 12 and Table 5, IC1 is maximal under El Niño conditions, positive NPGO (colder central North Pacific), slightly negative PDO, warmer (colder) mid-latitude (subtropical) Atlantic region, a dipolar structure in the South Pacific, and cold pools southward of 30°S in the Indian and Atlantic oceans. That is consistent with the oceanic influence of ENSO on decadal timescales [73]. Moreover, that agrees with the dynamical interactions between ENSO, PDO, and NPGO which have been discussed by several authors [74–76]. The quite relevant EOF11 pattern in IC1 is mostly correlated with the El Niño SST composite (asymmetric with respect to La Niña) at mid-latitudes in all oceans, during the December-February season [68, 77], which explains the above referred coherence of PC1 and PC11 during strong El Niños.The IC2 loading map (Figure12) mainly coincides with the symmetrical EOF6, that is, with the positive Atlantic Niño phase, while the remaining dominant EOF contributors to IC2, that is, EOF1, EOF7, and EOF9 (Table 5), practically cancel each other, especially in tropical Pacific El Niño region, despite the presence of non-null correlation in that region (Figure 11). Negative IC2 extremes (e.g., in 1912 and 1951) coincide with the strongest Atlantic Niñas in agreement with the positive skewness of PC6 (Table 2).The IC3 loading map pattern (Figure12) is quite similar to EOF5 (Central Pacific (CP) El Niño), combined with the NPGO (colder North Pacific for positive IC3). That is also visible in the correlation map (Figure 11). The negative skewness of IC3 (Table 5) is compatible with the primarily negative SST skewness in Central Pacific [78]. It supports the existence of strong Central Pacific La Niñas (strong negative values of IC3), as in 1943-1944 (see Figure 10), especially during the positive phase of AMO (see the negative loading factor of PC2 for IC3 in Table 5). The loading factors of PC4, PC5, and PC2 contributing to IC3 (Table 5) also reflect the tendency for CP-El Niño (CP-La Niña) during negative (positive) phases of AMO [79], which is visible through the symmetrical correlation signals of IC3 with SSTA (Figure 11), respectively, in CP and North Atlantic areas. The imprint of the characteristic South Pacific dipole shown in EOF10 is also visible in the loading map (Figure 12).The IC4 loading map (Figure12) is a mix of several EOFs. The inverse proportionality of IC4 to the PCs standard deviation makes EOFs 4, 5, 9, 10, and 11 the strongest contributors (Table 5 and Figure 2(a)) to the loading map of IC4 (Figure 12), showing that IC4 is highly positive when the Indian and Atlantic oceans are warm, while the eastern (western) North and South Pacific are warm (cold), including the tropical El Niño region. That is also reflected in the correlation map (Figure 11).The last IC5 loading map pattern (Figure12) is nearly symmetrical to EOF10 (Table 5), especially over the Austral oceans with a characteristic sequence of three dipoles, located over the Indian, Pacific, and Atlantic subtropical oceans. Negative extremes of IC5 (e.g., 1947, 1992/93) are also associated with a strong negative AMO and NPGO. The NE of IC5 result essentially from the highly positive kurtosis and long pdf tails of PC10 and PC2 (Table 2). Both the loading and correlation maps are quite similar.The results of ISA are summarized in Figure13. They suggest 5 ICs with the emergence of one dyadic link (IC3, IC5) with ITI = 0.015 and a triadic link (IC1, IC4, and IC5) with ITI = 0.024 and finally the strongest dyad (IC1, IC5) with an ITI = 0.045. This interaction is even larger than the self-NEs of ICs 3, 4, and 5. This interaction is mainly due to the nonlinear correlation: cor(IC1,IC52)=0.23, that is, extremes of IC5 occurring preferentially under PC1 > 0 (e.g., strong El Niños). A possible explanation is that PC4 (simulating NPGO) influence both IC1 and IC5. From the ensemble of interactions, only IC2 appears to be significantly independent of the remaining ICs. All the remaining ITIs remain at the noise level (Figure 13).Figure 13
As in Figure9 but for the ISA summary of SSTA data.The above analysis extracts highly negentropic ICs resulting primarily from independent combinations and coherency between PC extremes. An alternative ICA and ISA analysis is obtained with Gaussianized (followed by whitening) PCs. Further experiments (not shown) indicate that the total NE comes from cross cumulants only and ICs emerge exclusively from combinations of joint PC extremes. The two leading ICs are practically unchanged as compared to previous ICs; the third and fourth ICs are slightly mixed while the fifth IC is not significant, that is, remaining in the Gaussian manifold. In summary, a prior Gaussianization of PCs leads in general to a reduction, both of the dimension and of the strength of the non-Gaussian subspace.


### 3.5. Sensitivity to the PC Spanning Space

Here, we study the behavior of the non-Gaussian subspace and ICA/ISA outputs when we sequentially increase the metric space dimensionNmax using the leading whitened PCs. In the discussion below, it is worth mentioning that, given the direct sum: S+=S-⊕Sr between two orthogonal subspaces S-, Sr, the relationships between the associated non-Gaussian subspaces and their dimensions hold:(19)WngS+⊇WngS-⊕WngSr;NngS+≥NngS-+NngSr,which leads to the identity(20)WngS+·WngS-=NngS-.Hereafter, we focus on the case where S-, Sr are spanned, respectively, by the leading Nmax-N′ PCs and by the remaining N′≥1 PCs up to rank Nmax. In (19), the equality (19) holds when S-, Sr are statistically independent or equivalent when spaces S-, Sr have separable pdfs. Moreover, under statistical independence, the ISA/ICA may be done independently. However, the pdf separability is hardly verified on multiscale nonlinear natural systems as studied here which means that Wng(S+) intersects the Gaussian subspaces Wg(S-), Wg(Sr) of the parcel spaces while ICA/ISA is also changing when passing from S- to S+.An equivalent relationship to (19) for Wng-Ed, obtained with finite samples, is much harder to obtain since Wng-Ed may be undersized (oversized) by errors of type I (II) when considering the Gaussianity null hypothesis HG. This means that Wg-Ed(S+)·Wng-Ed(S-)≥0; that is, there may be some “leakage” of the “prior” non-Gaussian subspace (Wng-Ed(S-)), especially through the less non-Gaussian ICs, to the “posterior” Gaussian subspace Wg-Ed(S+).Next, we analyze how the ICA/ISA applied to SSTAs varies withNmax. To this end, we have considered the sequence Nmax=3,9,11,15,17 in which the number of ICs Nng-Ed, obtained by the method of Section 2.2.2 using a confidence level psig=90% (see, e.g., Figure 6), is Nng-Ed=1,4,5,6,7, respectively. The above Nmax are those for which Nng-Ed have shown a jump (of one or more dimensions) with respect to that obtained for Nmax-1, where more substantial changes are expected to occur in the ICA/ISA sources. For example, Nng-Ed=5 for 11≤Nmax≤14. Let us consider the corresponding sequence of embedding spaces: S-, S+, denoting by r1-,…,rN-- and r1+,…,rN++ the IC weighting vectors (sorted by decreasing self non-Gaussianity), spanning Wng-Ed(S-) and Wng-Ed(S+), respectively.To assess the impact on ICA/ISA, we have computed some inner-product-based diagnostics, plotted in Figure14. First, the measure P-+≡Wng-Ed(S+)·Wng-Ed(S-)=∑i=1N-∑j=1N+(ri-·rj+) while the difference P-+-N-≥0 yields how the non-Gaussian subspaces are preserved throughout the space sequence, that is, how closely (20) is preserved. Then, we identify which terms 0≤ri-·rj+≤1 are mostly contributing to that, typically above a threshold (e.g., 0.08). For dimensions within the range of the selected Nmax where Nng-Ed is kept constant, every ri- is very well projected onto a certain rj(i)+, leading to well tracked and robust ICs independently of the space dimension. For instance, in the sequence Nmax=11,12,13,14 (Nng-Ed=5), we get P-+=4.7,4.6,4.8, respectively, (quite close to the ideal value of 5) with very low leakage of non-Gaussianity. We also get the tracks of the five ICs: [IC1→IC10.96→IC10.91→IC10.94]; [IC2→IC30.99→IC30.98→IC30.95]; [IC3→IC20.83→IC20.77→IC20.99]; [IC4→IC40.97→IC40.98→IC40.97]; and [IC5→IC50.90→IC50.88→IC50.91], each one within brackets, where the subscript values stand for ri-·rj(i)+. We conclude that ICA and ISA are quite robust in the range Nmax = 11–14. Only the ranks of IC2 and IC3 were swapped in the step Nmax=11→12.Figure 14
Schematic of the ICs (rank and negentropy shown in rectangles), aligned vertically for eachNmax=3,9,11,15,17 (circles) with values of P-+ (bulk arrows) between the sequential PC spanning spaces. Values with the dominant ri-·rj+ are written next to the corresponding thin arrows. The red straight vertical lines indicate the significant dyadic links with corresponding NE interactivities.Now, in Figure14, we analyze what happens when a highly non-Gaussian or highly nonlinearly interacting PC comes into the analysis space, corresponding, for instance, to PC3, PC9, PC11, PC15, and PC17. EOF15 and EOF17 are both dominated by subtropical modes in the South Pacific and Indian oceans (not shown), nonlinearly interacting with precedent modes.In every sequenceS-→S+, the ICs may experience three situations: (1) the emergence of new ICs; (2) the “leakage” of non-Gaussianity to the forthcoming Gaussian subspace; and (3) the decoupling of sets of ICs (one or more) into sets of ICs of larger dimensions. In our case, the emergence of 3 new ICs appears when Nmax changes from 3 to 9 (see Figure 14). The “leakage” is more relevant with weaker ICs (low negentropy). The decoupling of ICs becomes frequent, mainly at higher dimensions, being the main mechanism of IC generation. For instance IC3 at Nmax=9 decouples into (IC4, IC5) with Nmax=11. Multiple interactivities, like the dyadic ones (expressed by red segments in Figure 14), may emerge on incoming ICs or transfer to the tracked ICs (e.g., the dyad (IC1, IC6) with Nmax=15 transfers to the dyad (IC1, IC7) with Nmax=17). There are some ICs that remain robust without any decoupling like the case of IC1 and IC3 with Nmax=11. In general, as Nmax increases, ICs project progressively onto low-ranked PCs dominated by smaller spatial scales and complex sharper features. In this situation, a hybrid criterion for the choice of the truncation Nmax (not explored here) may be favored by four factors: the negentropy and explained variance by ICs, the low dispersion of IC weighting loadings (either on PC or original basis), and a low Nng-Ed. That justifies the choice Nmax=11 as being the lowest dimension for which Nng-Ed=4.


## 3.1. Data Description and Processing

We consider the monthly SST data from the “Extended Reconstruction Sea Surface Temperature” (ERSST) version v3b product (https://www.ncdc.noaa.gov/oa/climate/research/sst/ersstv3.php), from the International Comprehensive Ocean–Atmosphere Data Set (ICOADS) release 2.4. The SST data are on 2°×2° 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 Nt=1224 monthly SST anomalies (SSTAs), on Np=8714 grid points. A EOF analysis is then applied [51] and we only use the space spanned by the leading Nmax=12 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 σTotal2=1709.4 (°C)2, the local variance at point P is σP2, with an average σLocal2≡σTotal2/[NPclat]~(0.50)2 (°C)2 (clat = average of cos⁡(latitude)) and we let the fractional local variance σ^P2≡σP2/σLocal2. 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 c^1 based on an AR(1) model.Figure 2
(a) % of explained variance confidence interval and (b) lag-one autocorrelation of the AR(1) fitting model of the leading 12 PCs.
(a)
(b)Figure3 shows the leading 12 EOFs which are described, with references in Table 1. Based on c^1, an estimate of the “average interval between independent realizations” is Δind=Nt/Ndof=fREd(Nt,c^1,N), (Figure 1(a)) which is not too different from (-1.)/ln⁡(c^1) (e.g., [53]).Table 1
EOF description with the main SSTA modes, along with corresponding available indices from literature and references.
EOF rank
Description with the main SSTA modes
Available index
PC/Index cor.
Refs.
Notes
1
El Niño (ENSO)
MEI
0.98
[56]
2
Atlantic Multi-Decadal Oscillation (AMO)
AMO
0.80
[57]
3
Pacific Decadal Oscillation (PDO): negative phase
PDO
-
0.56
[58, 59]
4
Northern Pacific Gyre Oscillation (NPGO)
NPGO
0.53
[60]
El Niño Modoki (EM): negative phase
EMI
-
0.42
[61]
Also known as Central Pacific or noncanonical El Niño
5
El Niño Modoki (EM): positive phase
EMI
0.57
[61]
Subtropical Indian Ocean Dipole: negative phase
SIOD
[62, 63]
First EOF mode of the SSTA in the Indian Ocean ranging in the 15°S–45°S band
6
Atlantic Niño: negative phase
ATL3
-
0.49
[64, 65]
South Atlantic Subtropical Dipole
SASD
[66]
First EOF coupled mode of the SSTA/wind in the South Atlantic
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
Zonally oriented dipole in South Pacific
[67, 68]
Second EOF of the South Pacific SSTA
7
Meridional sequence of dipoles crossing the North-South Pacific
South Pacific Ocean Dipole (SPOD)
SPOD
[67, 69]
First EOF of the South Pacific SSTA
8
North Atlantic tripole (NAT)
NAT
[70]
Tropical Atlantic SST dipole
[68]
9
South Pacific Ocean Dipole (SPOD)
SPOD
[67, 69]
First EOF of the South Pacific SSTA
10
Zonally oriented dipole in South Atlantic
[68]
Zonally oriented dipole in South Pacific
[67, 68]
Second EOF of the South Pacific SSTA
11
Meridional 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
12
Tripolar pattern in North Pacific
North Atlantic tripole (NAT)
NAT
[70]Table 2
% of explained variance, skewness, and kurtosis of the 12 leading PCs and the squared projections ofWng-Ed onto the individual whitened PCs for the analysis over N=11 PCs (see the end of Section 3.2).
PC1
PC2
PC3
PC4
PC5
PC6
PC7
PC8
PC9
PC10
PC11
PC12
%Var
17.6
5.3
4.6
3.8
3.3
2.7
2.5
2.2
2.2
2.0
1.9
1.8
Skew.
0.24
0.03
0.01
0.28
-
0.16
0.25
-
0.08
0.23
0.14
0.16
0.29
0.13
Kurt.
-
0.01
-
0.64
-
0.33
0.19
0.02
-
0.02
-
0.06
0.17
-
0.04
0.55
0.14
-
0.16
SqPro
0.65
0.53
0.16
0.68
0.50
0.56
0.29
0.05
0.37
0.81
0.39Figure 3
Maps of the 12 variance-leading, spatially normalized EOFs of the monthly SSTAs field, computed in the period 1910–2011. Ranked EOFs are read from left to right and top to bottom.Figure2(b) shows that Δind 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 c^1 values, that is, 0.86 (see Section 3.3).


## 3.2. Negentropy Distribution and Non-Gaussian Subspace Estimation

Non-Gaussianity and negentropy of a fieldx^ 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 y, through linear operators. Let us consider, for example, the local-basis coskewness tensor S(x^). 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 (P,Q,R):(17)Sx^F2≡∑P,Q,R=1NPSx^PQR2σ^P2σ^Q2σ^R2a^Pa^Qa^R=∑i,j,k=1NrSyijk2λ^iλ^jλ^k≃SyF2Nr3∑i,j,k=1Nrλ^iλ^jλ^k,where a^P stand for the relative area at grid point P and σ^P2 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 NEJ~Ed(y) can be decomposed using M~(y) (5). For each N=1,…,Nmax=12, we denote by y(N)=(Y1,…,YN) the set of N leading whitened PCs. Figure 4 shows the diagonal terms (M~kk[y(N)]), Figure 4(a) and the total NE J~Ed[y(N)]=Tr(M~[y(N)]), Figure 4(b). Note that, for a given N, the smaller values of M~kk[y(N)] are either the least non-Gaussian PCs or the less interactive PCs with the remaining PCs (i≠k), providing therefore good candidates for the Gaussian subspace of y(N). For example, with N=12, we have k=2,7,8,9, 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 M~22[y(N)] of the PC2 due to its small degrees of freedom.Figure 4
(a) NE components parcelM~kk[y(N)] (k in abscissas) for each N=1,…,12 (marked in grey shaded). (b) Total NE J~Ed[y(N)] for N=1,…,12. Values rejecting the null hypothesis of multi-Gaussianity at confidence levels <70%, 70%–90%, 90%–95%, and >95% are marked, respectively, by “+,” open circle, small solid circle, and large solid circle.
(a)
(b)BiasesE{J~Ed[y(N)]}-JEd[y(N)] may be approximated as E(J~Ed:G)=fREd(Nt,c^1,N)PJEd(N)/Nt (see (11)) with c^1≈0.86, leading to fREd(Nt=1224,c^1=0.86)~3.3 (see Figure 1(a)). For instance, for N=12, J~Ed=3.1, 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.Figure4(a) shows that, for all N, the largest NE component is M~11[y(N)], 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:M~kk[y(N)]-M~kk[y(N-1)] (Figure 4(a)), one may detect which interactivities mostly contribute to J~Ed[y(N)]. The largest differences are for (k,N) pairs: {(1,3),(1,5),(1,10),(1,11), (3,11),(4,10),(4,11),(5,10),(5,11),(6,11)} in consistency with the high confidence level of M~kk[y(11)] for k=1,3–6,10,11. Those jumps are often due to particular joint cumulants in the 2D subspaces spanned by (PCk,PCN), responsible for particular types of joint non-Gaussianity, for example, scatter-plots resulting for instance from the occurrence of joint extremes. For instance, Yk2Yi¯>0, is favored for the sign combinations (-,+) and (+,+) of (Yk,Yi). As an illustration, Figures 5(a) and 5(b) show, respectively, scatter-plots of (Y1g,Y3g) and (Y1g,Y11g) (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 3 of [50]) during which a positive PC11 and a negative PDO phase have occurred.Figure 5
Scatter-plots of PC1 versus PC3 (a) and PC1 versus PC11 (b). The marginal pdfs have been Gaussianized.
(a)
(b)SinceM~12,12[y(12)] 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 Hg test to limit N [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 vectory(11) and the associated vector z=Z1,…,Z11′, yielding J~Ed[z]=2.55, and a bias of 1.73, leaving ~0.82 for “true non-Gaussianity.” Figure 6 shows the NE components (or SVs) M~kk(z) (k=1,…,N), along with the threshold quantiles 90% and 95% of Qpsig[M~kk(bk)] obtained from an ensemble of 1000 Gaussian surrogates. For comparison, we also show the SVs M~kk(zs) (open squares) of surrogate realizations zs of z. The SVs of zs all lie below the above thresholds, supporting the consistency of the Gaussianity test. Figure 6 also suggests that the estimated Gaussian subspace Wg-Ed of the SSTA is spanned by the last 6 SVecs (with a confidence level psig=95%) and so Ng-Ed=6 and, Nng-Ed=5 for the non-Gaussian subspace Wng-Ed. The total NE restricted to zng≡Z1,…,Z5′ is J~Ed[zng]=1.11 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 M~kk(zng) (open black circles), with their low values compared to M~kk(z), k=1,…,5, are also shown in Figure 6.Figure 6
Ranked SVsM~kk(z), k=1,…,N=11, contributing to NE, both from real SSTA data (solid squares) and from Gaussian surrogate realizations (open squares) and 90% (red), 95% (blue) quantiles of Qpsig[M~kk(bk)] obtained from Gaussian surrogates after rank k, to be compared with M~kk(z). The trailing 6 SVs is not significant at psig=95%, leading to Nng-Ed=5. The NE components restricted to the non-Gaussian subspace are marked as open circles.The squared projections ofWng-Ed onto each of the 11 whitened PCs are given by the inner product (Appendix B) Wng-Ed·ei=∑k=15Ri,k2, for i=1,…,N=11, (Table 2), where ei is the (unitary) ith vector. By definition, those values add up to Nng-Ed=5. 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 ∑i=111λ^i∑k=15Ri,k2 = 24.7% and 23.3%, for the non-Gaussian and Gaussian subspaces.


## 3.3. Application to Synthetic Non-Gaussian Data

We generateNmax=12 AR(1) independent series Gk(t), t=1,…,Nt, k=1,…,Nmax. The whitened and uncorrelated components of the non-Gaussian vector y′=(Y1,…,YN) are computed following(18)Yk=Gk;k=1,…,N≤Nmax,k≠2,5Y2=σ2-1/2G12-1+1-σ21/2G2;Y5=σG3G4+1-σ21/2G5with σ2∈[0,1] measuring the “strength” of non-Gaussianity. So y is composed of Ng=N-5 Gaussian scalar sources spanning Wg and Nng=5 non-Gaussian variables spanning Wng. The latter is composed of one dyad y′(1)=(Y1,Y2) and a non-Gaussian triad y′(2)=(Y3,Y4,Y5).The next step consists of testing the null hypothesisHG of global Gaussianity (i.e., σ2=0) and evaluating J~Ed(y) by varying σ2 in [0,1] and Ng in {0,…,7}.Figure7(a) shows the Monte-Carlo average E[J~Ed(y)] obtained from 1000 surrogates. As expected, the mean NE increases with σ2 (nonlinearity effect) and Ng (dimensionality effect), though at a much smaller rate. The bias E[J~Ed:G(y)] is the value reached at σ2=0 in agreement with (11).Figure 7
(a) Monte-Carlo average of the sample negentropy estimatorJ~Ed(y) as a function of σ2 and Ng. (b) Contour plot of the fraction F of times for which the null hypothesis HG is rejected at psig=99% (1% significance level) as a function of σ2and Ng. The contour F=0.99 can mark σmin2 at the 99% confidence level.
(a)
(b)


### 3.3.1. Gaussian Subspace Recovery

We investigate the strength ofJ~Ed(y) as a function of σ2 and Ng. The null hypothesis HG is rejected at the confidence level psig if J~Ed(y)>[PJEd(N)Nt-1]fREd-Psig using fREd-Psig (Figure 1(b)). The Monte-Carlo fraction F of HG rejection is shown in Figure 7(b) for psig=99%. The contour line at F=0.99 provides the graph σmin2(Ng) of σ2 for the recovery of the non-Gaussian sources. As expected, σmin2 increases with Ng, starting at σmin2=20% and at Ng=0 and reaching σmin2~45% at Ng=6. This means that the larger N, the stronger (through σ2) the non-Gaussian embedded sources. Moreover, σmin2 diminishes under a less conservative test (lower psig) 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., psig=90%) [47].


### 3.3.2. ICA and ISA

Here, we apply ICA and ISA forN=11 (Ng=6) and compute the skill of source separation in two cases: (a) σ2=20%, high noise, and (b) σ2=40%, intermediate noise, yielding, respectively, total NE of 1.79 and 2.57.Figures8(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: Nng-Ed=5 with psig=90% 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 Wng-Ed·Wng (compared to 5 in the ideal recovery). Those products are, respectively, 3.99 and 4.49 for cases (a) and (b). Lower values of σ2 lead to a subestimation of Nng and to a worse estimation of Wng.Figure 8
Same as Figure6 but for the for synthetic data with σ2=20% (a) and σ2=40% (b).
(a)
(b)Now, by applying ICA (see Section2.2.3 and Appendix A) we get ICs making the vector zICA=(RICA′R′)y=Z1,ICA,…,ZNNg-Ed,ICA′. Tables 3 and 4 give, respectively, for cases (a) and (b), a summary of ICA with the sorted maximized self-NEs J~Ed(Zi,ICA), the NE components M~ii(zICA), and the IC vector loadings in RRICA. 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 zEd-dyad is spanned by IC1 and IC4 whereas the estimated triad zEd-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 σ2. 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).Table 3
Contributions to the ICA contrast function (2nd column). NE components (3rd column) and IC vector loadings on the basis of original variables (18) for the high noise case. Components in which ICs are most projected are in bold font.
IC rank
J
~
E
d
(
Z
i
,
I
C
A
)
M
~
i
i
(
z
I
C
A
)
Y
1 (dyad)
Y
2 (dyad)
Y
3 (triad)
Y
4 (triad)
Y
5 (triad)
IC1:i=1
0.127
0.162
0.76
-
0.58
−0.02
0.01
0.16
IC2:i=2
0.066
0.125
−0.07
0.07
0.47
-
0.56
0.56
IC3:i=3
0.042
0.102
−0.01
−0.08
-
0.38
-
0.68
-
0.42
IC4:i=4
0.014
0.053
-
0.41
-
0.52
0.21
0.29
0.01
IC5:i=5
0.011
0.063
−0.03
0.02
-
0.68
0.07
0.57
Total
0.260
0.506Table 4
Table3 but for the case of intermediate noise (b).
IC rank
J
~
E
d
(
Z
i
,
I
C
A
)
M
~
i
i
(
z
I
C
A
)
Y
1 (dyad)
Y
2 (dyad)
Y
3 (triad)
Y
4 (triad)
Y
5 (triad)
IC1:i=1
0.410
0.454
0.70
-
0.70
−0.03
0.00
0.08
IC2:i=2
0.256
0.372
0.00
0.05
0.45
0.65
0.58
IC3:i=3
0.109
0.257
0.03
−0.03
-
0.30
0.72
-
0.55
IC4:i=4
0.033
0.078
0.59
0.59
−0.17
−0.07
0.00
IC5:i=5
0.031
0.114
−0.04
0.04
-
0.76
0.11
0.54
Total
0.840
1.277The final step consists of merging the ICs as outlined in Section2.2.3. Figure 9 presents the results for both cases, where the self-interactivities (16) are the most relevant ITIs due to their maximization by ICA (except IC5 in case (a)). The cumulated NE reaches the unbiased NE value (dashed line) near the threshold of the ITI acceptance (scaled-ITI > 1). In case (a) (Figure 9(a)), the accepted ITIs are for ICs 1–5, and dyadic links (IC3, IC2), (IC5, IC3) yielding the emergence of the true triad (IC2, IC3, and IC5) (see Table 3). The marginally accepted ITIs holds for (IC4, IC3, and IC2) (false triad) and (IC1, IC4) (true dyad). Excluded ITIs remain at the level of the NE bias. For smaller values of σ2 (not shown), the recovered sources may be mixed, partly restored, or even not detected.Figure 9
Sorted scaled-interactivities by decreasing order (solid circles), ITIs (solid squares) with cumulated ITIs (open crossed circles) and unbiased total NE (horizontal dashed line). The sorted IC groups are labeled (e.g., 32 means group of IC2 and IC3). The 95% confidence levels of the selected ITIs are those for which the cumulated ITIs are below the dashed line. (a) Case (a): high noise and (b) case (b): intermediate noise.
(a)
(b)For case (b) (Figure9(b)) the source recovery improves. In fact, the significant accepted ITIs are for ICs 1–5, the links for (IC5, IC3), (IC5, IC2), (IC3, IC2) with the true emerging triad (IC2, IC3, and IC5), and for the true dyad (IC4, IC1). The sum of ITIs associated with ICA and ISA is, respectively, 0.84 and 1.13, showing clearly the advantage of ISA in this case.


## 3.3.1. Gaussian Subspace Recovery

We investigate the strength ofJ~Ed(y) as a function of σ2 and Ng. The null hypothesis HG is rejected at the confidence level psig if J~Ed(y)>[PJEd(N)Nt-1]fREd-Psig using fREd-Psig (Figure 1(b)). The Monte-Carlo fraction F of HG rejection is shown in Figure 7(b) for psig=99%. The contour line at F=0.99 provides the graph σmin2(Ng) of σ2 for the recovery of the non-Gaussian sources. As expected, σmin2 increases with Ng, starting at σmin2=20% and at Ng=0 and reaching σmin2~45% at Ng=6. This means that the larger N, the stronger (through σ2) the non-Gaussian embedded sources. Moreover, σmin2 diminishes under a less conservative test (lower psig) 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., psig=90%) [47].


## 3.3.2. ICA and ISA

Here, we apply ICA and ISA forN=11 (Ng=6) and compute the skill of source separation in two cases: (a) σ2=20%, high noise, and (b) σ2=40%, intermediate noise, yielding, respectively, total NE of 1.79 and 2.57.Figures8(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: Nng-Ed=5 with psig=90% 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 Wng-Ed·Wng (compared to 5 in the ideal recovery). Those products are, respectively, 3.99 and 4.49 for cases (a) and (b). Lower values of σ2 lead to a subestimation of Nng and to a worse estimation of Wng.Figure 8
Same as Figure6 but for the for synthetic data with σ2=20% (a) and σ2=40% (b).
(a)
(b)Now, by applying ICA (see Section2.2.3 and Appendix A) we get ICs making the vector zICA=(RICA′R′)y=Z1,ICA,…,ZNNg-Ed,ICA′. Tables 3 and 4 give, respectively, for cases (a) and (b), a summary of ICA with the sorted maximized self-NEs J~Ed(Zi,ICA), the NE components M~ii(zICA), and the IC vector loadings in RRICA. 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 zEd-dyad is spanned by IC1 and IC4 whereas the estimated triad zEd-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 σ2. 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).Table 3
Contributions to the ICA contrast function (2nd column). NE components (3rd column) and IC vector loadings on the basis of original variables (18) for the high noise case. Components in which ICs are most projected are in bold font.
IC rank
J
~
E
d
(
Z
i
,
I
C
A
)
M
~
i
i
(
z
I
C
A
)
Y
1 (dyad)
Y
2 (dyad)
Y
3 (triad)
Y
4 (triad)
Y
5 (triad)
IC1:i=1
0.127
0.162
0.76
-
0.58
−0.02
0.01
0.16
IC2:i=2
0.066
0.125
−0.07
0.07
0.47
-
0.56
0.56
IC3:i=3
0.042
0.102
−0.01
−0.08
-
0.38
-
0.68
-
0.42
IC4:i=4
0.014
0.053
-
0.41
-
0.52
0.21
0.29
0.01
IC5:i=5
0.011
0.063
−0.03
0.02
-
0.68
0.07
0.57
Total
0.260
0.506Table 4
Table3 but for the case of intermediate noise (b).
IC rank
J
~
E
d
(
Z
i
,
I
C
A
)
M
~
i
i
(
z
I
C
A
)
Y
1 (dyad)
Y
2 (dyad)
Y
3 (triad)
Y
4 (triad)
Y
5 (triad)
IC1:i=1
0.410
0.454
0.70
-
0.70
−0.03
0.00
0.08
IC2:i=2
0.256
0.372
0.00
0.05
0.45
0.65
0.58
IC3:i=3
0.109
0.257
0.03
−0.03
-
0.30
0.72
-
0.55
IC4:i=4
0.033
0.078
0.59
0.59
−0.17
−0.07
0.00
IC5:i=5
0.031
0.114
−0.04
0.04
-
0.76
0.11
0.54
Total
0.840
1.277The final step consists of merging the ICs as outlined in Section2.2.3. Figure 9 presents the results for both cases, where the self-interactivities (16) are the most relevant ITIs due to their maximization by ICA (except IC5 in case (a)). The cumulated NE reaches the unbiased NE value (dashed line) near the threshold of the ITI acceptance (scaled-ITI > 1). In case (a) (Figure 9(a)), the accepted ITIs are for ICs 1–5, and dyadic links (IC3, IC2), (IC5, IC3) yielding the emergence of the true triad (IC2, IC3, and IC5) (see Table 3). The marginally accepted ITIs holds for (IC4, IC3, and IC2) (false triad) and (IC1, IC4) (true dyad). Excluded ITIs remain at the level of the NE bias. For smaller values of σ2 (not shown), the recovered sources may be mixed, partly restored, or even not detected.Figure 9
Sorted scaled-interactivities by decreasing order (solid circles), ITIs (solid squares) with cumulated ITIs (open crossed circles) and unbiased total NE (horizontal dashed line). The sorted IC groups are labeled (e.g., 32 means group of IC2 and IC3). The 95% confidence levels of the selected ITIs are those for which the cumulated ITIs are below the dashed line. (a) Case (a): high noise and (b) case (b): intermediate noise.
(a)
(b)For case (b) (Figure9(b)) the source recovery improves. In fact, the significant accepted ITIs are for ICs 1–5, the links for (IC5, IC3), (IC5, IC2), (IC3, IC2) with the true emerging triad (IC2, IC3, and IC5), and for the true dyad (IC4, IC1). The sum of ITIs associated with ICA and ISA is, respectively, 0.84 and 1.13, showing clearly the advantage of ISA in this case.


## 3.4. Application of ICA and ISA to SSTA Data

Using the estimated non-Gaussian subspace (Section3.2), we apply here ICA to the leading NICA=5 SVecs of the SSTA dataset (see Figure 6). Table 5 gives a summary of the application with the characteristics of the ranked ICs Zi,ICA, i=1,…,5, by decreasing self-NEs. They explain, respectively, 7.1%, 5.5%, 4.3%, 4.1%, and 3.7% of the total variance.Table 5
Sorted ICsi=1,…,5, for SSTA data showing self-NEs contributing to the ICA contrast function (2nd column), NE components (3rd column), and loading factors (RRICA)k,i (i=1,…,5), of the leading ICs on the PC basis. Loading factors, sorted by absolute value, explaining up to 85% of the unitary norm are in bold font. The last two columns give the skewness and kurtosis of ICs.
IC rank
J
~
Ed
(
Z
i
,
I
C
A
)
M
~
i
i
(
z
I
C
A
)
PC1
PC2
PC3
PC4
PC5
PC6
PC7
PC8
PC9
PC10
PC11
Sk.
Kurt.
IC1:i=1
0.694
0.747
0.53
0.05
0.34
0.40
−0.17
0.24
0.06
0.12
−0.10
0.18
0.52
1.48
4.95
IC2:i=2
0.099
0.137
0.44
0.13
−0.14
0.02
0.03
-
0.63
-
0.35
−0.10
-
0.44
0.13
−0.15
−0.62
1.79
IC3:i=3
0.042
0.091
0.24
-
0.44
0.11
0.35
0.55
-
0.27
0.24
0.03
0.28
-
0.29
−0.12
−0.44
1.12
IC4:i=4
0.042
0.074
−0.26
0.35
−0.07
0.47
0.36
0.12
0.28
0.10
-
0.29
0.42
-
0.27
−0.24
1.32
IC5:i=5
0.016
0.064
−0.21
0.44
0.06
0.42
−0.18
−0.12
−0.14
0.10
−0.04
-
0.70
0.03
−0.32
0.58
Total
0.893
1.113Most of the NE components come from the self-NEs (an effect of ICA), through their high skewness and positive kurtosis or leptokurtosis (Table5), leading to fat long pdf tails. Extreme values of about 3–5 standard deviations are often reached as seen by the IC time-series (Figure 10). Under those conditions, ICs are basically aggregative and cooperative combinations of single and joint PC extremes. Quasi-independence and noncorrelation between ICs means that occurrence of extremes among different ICs tends to be asynchronous (see Figure 10) and therefore extreme phenomena are independent.Figure 10
IC time-series for the period 1910–2011.Figure11 shows the correlation maps between the ICs and SSTAs. An alternative way to interpret ICs is through the patterns displayed in the corresponding normalized loading maps (Figure 12) obtained from combinations of EOFs (Figure 3). Correlation and loading map values at the grid point P for the ith IC are weighted by the PCs standard deviation and its inverse, respectively, ∑k=1N=11Dk,Pek(P)(RRICA)k,i where Dk,P is λ^k1/2/σ^P (Figure 11) and λ^k-1/2 (Figure 12), respectively, and ek(P) is the kth EOF at P (Figure 3). By definition, ICs are spatial inner products between correspondent normalized loading maps (Figure 12) and the SSTA field. Consequently, maximal positive (negative) ICs on a given month are favored by highly positive (negative) pattern correlations between their loading maps and the SSTA field. Note, however, that the patterns shown in Figures 11 and 12 are not orthogonal. On loading maps, the weighting gives rise to the impact on ICs of PCs of small explained variance. Because of this weighting, the correlation maps (Figure 11) often display much broader and large-scale patterns than loading maps (Figure 12) [15], though both share the same main features when the range of PC variances contributing to a given IC is not too wide, which is the case here.Figure 11
Correlation maps between ICs and SSTA field.Figure 12
Spatially normalized loading maps of the 5 sorted ICs.The leading IC1 explains about 2/3 of the total NE with the remaining 4 being almost evenly distributed (Table5). Further experiments show that IC1 is quite robust with the use of NICA<5 leading SVecs with the remaining ICs being slight mixtures of those obtained with NICA=5. The above mixture of PCs suggests the existence of a global nonseparable source thanks to the nonlinear interaction across scales. This strongly suggests that ICs, obtained as blind source separation, are convenient tools for identifying nonlinear variability. Dynamical constraints may be included explicitly, for example, via nonlinear dynamical modes [71, 72].The leading IC1 reaches values of +5 standard deviations in the years 1983 and 1997 of strong El Niños, during which PC1, PC3, PC4, and PC11 are synchronised (see also Figures5(a) and 5(b) and Table 5). The coherence between those four PCs may sometimes be partial.The correlation pattern for IC1 (Figure11) has an El Niño pattern shape due to the high PC1 amplitude (Table 5) and PC1 variance whereas the loading map (Figure 12) is much correlated with EOF11, showing more local features. Therefore, from Figure 12 and Table 5, IC1 is maximal under El Niño conditions, positive NPGO (colder central North Pacific), slightly negative PDO, warmer (colder) mid-latitude (subtropical) Atlantic region, a dipolar structure in the South Pacific, and cold pools southward of 30°S in the Indian and Atlantic oceans. That is consistent with the oceanic influence of ENSO on decadal timescales [73]. Moreover, that agrees with the dynamical interactions between ENSO, PDO, and NPGO which have been discussed by several authors [74–76]. The quite relevant EOF11 pattern in IC1 is mostly correlated with the El Niño SST composite (asymmetric with respect to La Niña) at mid-latitudes in all oceans, during the December-February season [68, 77], which explains the above referred coherence of PC1 and PC11 during strong El Niños.The IC2 loading map (Figure12) mainly coincides with the symmetrical EOF6, that is, with the positive Atlantic Niño phase, while the remaining dominant EOF contributors to IC2, that is, EOF1, EOF7, and EOF9 (Table 5), practically cancel each other, especially in tropical Pacific El Niño region, despite the presence of non-null correlation in that region (Figure 11). Negative IC2 extremes (e.g., in 1912 and 1951) coincide with the strongest Atlantic Niñas in agreement with the positive skewness of PC6 (Table 2).The IC3 loading map pattern (Figure12) is quite similar to EOF5 (Central Pacific (CP) El Niño), combined with the NPGO (colder North Pacific for positive IC3). That is also visible in the correlation map (Figure 11). The negative skewness of IC3 (Table 5) is compatible with the primarily negative SST skewness in Central Pacific [78]. It supports the existence of strong Central Pacific La Niñas (strong negative values of IC3), as in 1943-1944 (see Figure 10), especially during the positive phase of AMO (see the negative loading factor of PC2 for IC3 in Table 5). The loading factors of PC4, PC5, and PC2 contributing to IC3 (Table 5) also reflect the tendency for CP-El Niño (CP-La Niña) during negative (positive) phases of AMO [79], which is visible through the symmetrical correlation signals of IC3 with SSTA (Figure 11), respectively, in CP and North Atlantic areas. The imprint of the characteristic South Pacific dipole shown in EOF10 is also visible in the loading map (Figure 12).The IC4 loading map (Figure12) is a mix of several EOFs. The inverse proportionality of IC4 to the PCs standard deviation makes EOFs 4, 5, 9, 10, and 11 the strongest contributors (Table 5 and Figure 2(a)) to the loading map of IC4 (Figure 12), showing that IC4 is highly positive when the Indian and Atlantic oceans are warm, while the eastern (western) North and South Pacific are warm (cold), including the tropical El Niño region. That is also reflected in the correlation map (Figure 11).The last IC5 loading map pattern (Figure12) is nearly symmetrical to EOF10 (Table 5), especially over the Austral oceans with a characteristic sequence of three dipoles, located over the Indian, Pacific, and Atlantic subtropical oceans. Negative extremes of IC5 (e.g., 1947, 1992/93) are also associated with a strong negative AMO and NPGO. The NE of IC5 result essentially from the highly positive kurtosis and long pdf tails of PC10 and PC2 (Table 2). Both the loading and correlation maps are quite similar.The results of ISA are summarized in Figure13. They suggest 5 ICs with the emergence of one dyadic link (IC3, IC5) with ITI = 0.015 and a triadic link (IC1, IC4, and IC5) with ITI = 0.024 and finally the strongest dyad (IC1, IC5) with an ITI = 0.045. This interaction is even larger than the self-NEs of ICs 3, 4, and 5. This interaction is mainly due to the nonlinear correlation: cor(IC1,IC52)=0.23, that is, extremes of IC5 occurring preferentially under PC1 > 0 (e.g., strong El Niños). A possible explanation is that PC4 (simulating NPGO) influence both IC1 and IC5. From the ensemble of interactions, only IC2 appears to be significantly independent of the remaining ICs. All the remaining ITIs remain at the noise level (Figure 13).Figure 13
As in Figure9 but for the ISA summary of SSTA data.The above analysis extracts highly negentropic ICs resulting primarily from independent combinations and coherency between PC extremes. An alternative ICA and ISA analysis is obtained with Gaussianized (followed by whitening) PCs. Further experiments (not shown) indicate that the total NE comes from cross cumulants only and ICs emerge exclusively from combinations of joint PC extremes. The two leading ICs are practically unchanged as compared to previous ICs; the third and fourth ICs are slightly mixed while the fifth IC is not significant, that is, remaining in the Gaussian manifold. In summary, a prior Gaussianization of PCs leads in general to a reduction, both of the dimension and of the strength of the non-Gaussian subspace.


## 3.5. Sensitivity to the PC Spanning Space

Here, we study the behavior of the non-Gaussian subspace and ICA/ISA outputs when we sequentially increase the metric space dimensionNmax using the leading whitened PCs. In the discussion below, it is worth mentioning that, given the direct sum: S+=S-⊕Sr between two orthogonal subspaces S-, Sr, the relationships between the associated non-Gaussian subspaces and their dimensions hold:(19)WngS+⊇WngS-⊕WngSr;NngS+≥NngS-+NngSr,which leads to the identity(20)WngS+·WngS-=NngS-.Hereafter, we focus on the case where S-, Sr are spanned, respectively, by the leading Nmax-N′ PCs and by the remaining N′≥1 PCs up to rank Nmax. In (19), the equality (19) holds when S-, Sr are statistically independent or equivalent when spaces S-, Sr have separable pdfs. Moreover, under statistical independence, the ISA/ICA may be done independently. However, the pdf separability is hardly verified on multiscale nonlinear natural systems as studied here which means that Wng(S+) intersects the Gaussian subspaces Wg(S-), Wg(Sr) of the parcel spaces while ICA/ISA is also changing when passing from S- to S+.An equivalent relationship to (19) for Wng-Ed, obtained with finite samples, is much harder to obtain since Wng-Ed may be undersized (oversized) by errors of type I (II) when considering the Gaussianity null hypothesis HG. This means that Wg-Ed(S+)·Wng-Ed(S-)≥0; that is, there may be some “leakage” of the “prior” non-Gaussian subspace (Wng-Ed(S-)), especially through the less non-Gaussian ICs, to the “posterior” Gaussian subspace Wg-Ed(S+).Next, we analyze how the ICA/ISA applied to SSTAs varies withNmax. To this end, we have considered the sequence Nmax=3,9,11,15,17 in which the number of ICs Nng-Ed, obtained by the method of Section 2.2.2 using a confidence level psig=90% (see, e.g., Figure 6), is Nng-Ed=1,4,5,6,7, respectively. The above Nmax are those for which Nng-Ed have shown a jump (of one or more dimensions) with respect to that obtained for Nmax-1, where more substantial changes are expected to occur in the ICA/ISA sources. For example, Nng-Ed=5 for 11≤Nmax≤14. Let us consider the corresponding sequence of embedding spaces: S-, S+, denoting by r1-,…,rN-- and r1+,…,rN++ the IC weighting vectors (sorted by decreasing self non-Gaussianity), spanning Wng-Ed(S-) and Wng-Ed(S+), respectively.To assess the impact on ICA/ISA, we have computed some inner-product-based diagnostics, plotted in Figure14. First, the measure P-+≡Wng-Ed(S+)·Wng-Ed(S-)=∑i=1N-∑j=1N+(ri-·rj+) while the difference P-+-N-≥0 yields how the non-Gaussian subspaces are preserved throughout the space sequence, that is, how closely (20) is preserved. Then, we identify which terms 0≤ri-·rj+≤1 are mostly contributing to that, typically above a threshold (e.g., 0.08). For dimensions within the range of the selected Nmax where Nng-Ed is kept constant, every ri- is very well projected onto a certain rj(i)+, leading to well tracked and robust ICs independently of the space dimension. For instance, in the sequence Nmax=11,12,13,14 (Nng-Ed=5), we get P-+=4.7,4.6,4.8, respectively, (quite close to the ideal value of 5) with very low leakage of non-Gaussianity. We also get the tracks of the five ICs: [IC1→IC10.96→IC10.91→IC10.94]; [IC2→IC30.99→IC30.98→IC30.95]; [IC3→IC20.83→IC20.77→IC20.99]; [IC4→IC40.97→IC40.98→IC40.97]; and [IC5→IC50.90→IC50.88→IC50.91], each one within brackets, where the subscript values stand for ri-·rj(i)+. We conclude that ICA and ISA are quite robust in the range Nmax = 11–14. Only the ranks of IC2 and IC3 were swapped in the step Nmax=11→12.Figure 14
Schematic of the ICs (rank and negentropy shown in rectangles), aligned vertically for eachNmax=3,9,11,15,17 (circles) with values of P-+ (bulk arrows) between the sequential PC spanning spaces. Values with the dominant ri-·rj+ are written next to the corresponding thin arrows. The red straight vertical lines indicate the significant dyadic links with corresponding NE interactivities.Now, in Figure14, we analyze what happens when a highly non-Gaussian or highly nonlinearly interacting PC comes into the analysis space, corresponding, for instance, to PC3, PC9, PC11, PC15, and PC17. EOF15 and EOF17 are both dominated by subtropical modes in the South Pacific and Indian oceans (not shown), nonlinearly interacting with precedent modes.In every sequenceS-→S+, the ICs may experience three situations: (1) the emergence of new ICs; (2) the “leakage” of non-Gaussianity to the forthcoming Gaussian subspace; and (3) the decoupling of sets of ICs (one or more) into sets of ICs of larger dimensions. In our case, the emergence of 3 new ICs appears when Nmax changes from 3 to 9 (see Figure 14). The “leakage” is more relevant with weaker ICs (low negentropy). The decoupling of ICs becomes frequent, mainly at higher dimensions, being the main mechanism of IC generation. For instance IC3 at Nmax=9 decouples into (IC4, IC5) with Nmax=11. Multiple interactivities, like the dyadic ones (expressed by red segments in Figure 14), may emerge on incoming ICs or transfer to the tracked ICs (e.g., the dyad (IC1, IC6) with Nmax=15 transfers to the dyad (IC1, IC7) with Nmax=17). There are some ICs that remain robust without any decoupling like the case of IC1 and IC3 with Nmax=11. In general, as Nmax increases, ICs project progressively onto low-ranked PCs dominated by smaller spatial scales and complex sharper features. In this situation, a hybrid criterion for the choice of the truncation Nmax (not explored here) may be favored by four factors: the negentropy and explained variance by ICs, the low dispersion of IC weighting loadings (either on PC or original basis), and a low Nng-Ed. That justifies the choice Nmax=11 as being the lowest dimension for which Nng-Ed=4.


## 4. Summary and Conclusion

We present in this paper an approximation of the multivariate negentropy (NE) motivated by the Edgeworth expansion of the joint probability distribution, using a linear combination of Frobenius norms of the joint cumulant tensors of orderp≥3. This form is especially sensitive to the extent and frequency of extremes. The NE approximation is based on the trace of a positive definite matrix M resulting from the high-order singular value decomposition of an appropriate contracted form of the coskewness and cokurtosis joint tensors in consistency with the exact negentropy. First, it is a tensor invariant, that is, invariant under orthogonal rotations. Second, when variables are whitened, a version of the “Separation Source Theorem” holds; that is, the approximated NE may be explicitly decomposed into positive terms measuring self-NEs (due to skewness/kurtosis of marginals) and NEs resulting from nonlinear associations between variables. This permits the extraction of statistically independent components (ICs) from observations through the maximization of the sum of marginal NEs. The decomposition is further generalized to multivariate sources under the Independent Subspace Analysis (ISA), through the formation of groups of disjoint ICs, associated with positive NE contributors (or interactivities). The jointly Gaussian ICs span the Gaussian subspace, being simply approximated by the kernel subspace of M.When cumulants are estimated from finite or serially correlated samples, the NE estimator bias is positive and grows with the dimension (Curse of Dimensionality). Therefore, the estimated NE of a Gaussian distributed finite multivariate sample is prone to produce “false non-Gaussianity,” hence the importance of Gaussian subspace identification. This subspace is estimated from the trailing singular vectors ofM for which the null hypothesis of Gaussianity cannot be rejected. The search of non-Gaussian ICs or multiple sources is sought in the non-Gaussian subspace.The method is then applied to real data consisting of 11 leading PCs obtained from 102 years of the SST anomaly (SSTA) field over the global ocean, equatorward of 65°. To test the effectiveness of the method, prior to using the SSTA PCs, we have used synthetic data with a prescribed amount non-Gaussianity and sharing the same dimensions, sample size, lagged-one autocorrelation, and the level of negentropy as the SSTA PCs. In order to assure the recovery of the “true non-Gaussianity,” the experiment suggests a minimum of variance of the true non-Gaussian sources, which increases as the sample gets smaller and the dimension of the embedding space gets larger.Application of the method to SSTA data based on the leading 12 PCs (the last one going practically to the Gaussian subspace), explaining more than 50% of the monthly variance, yields 5 ICs. The most prominent or most negentropic IC (or IC1) represents El Niño conditions combined with negative PDO, positive NPGO, and composite impacts of El Niño (asymmetric with respect to La Niña) on oceanic mid-latitudes, poleward of 30°. The second component, IC2, reveals a loading map mostly projecting onto the Atlantic Niño and the AMO pattern. The third component, IC3, reflects the tendency of Central Pacific El Niños during the negative phase of AMO. The fourth component, IC4, combines several EOFs and the last component, IC5, is dominated by EOF10, being positive when the SSTA field in the South Pacific exhibits a zonally oriented dipole. All ICs are basically independent combinations of PC extremes and suggest independent responses of the SST field to different forcings. Only IC2 appears to be significantly isolated from the other components. The remaining components have some statistical links dominated primarily by the dyad (IC1, IC5), possibly associated with the sharing of the NPGO effect, via PC4, among those ICs.The identification of possible independent physical forcings, responsible for the independent sources may bring some benefits, concerning, namely, the forecasting and modelling of the SSTA variability and its impacts and teleconnections on relevant surface variables (e.g., precipitation, temperature, moisture, and wind). That opens the possibility of building phenomenological, dynamical, and stochastic models for the different ICs and independent subspaces.



---

*Source: 3076810-2017-08-22.xml*