Complexity

Volume 2017 (2017), Article ID 3076810, 23 pages

https://doi.org/10.1155/2017/3076810

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

^{1}Instituto Dom Luiz (IDL), Faculdade de Ciências, Universidade de Lisboa, 1749-016 Lisbon, Portugal^{2}Department of Meteorology, Stockholm University, Stockholm, Sweden

Correspondence should be addressed to Carlos A. L. Pires

Received 23 May 2017; Accepted 10 July 2017; Published 22 August 2017

Academic Editor: Davide Faranda

Copyright © 2017 Carlos A. L. Pires and Abdel Hannachi. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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.

#### 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. Section 2 presents the separation method and corresponding statistical tests. Section 3 presents application to the SST anomaly data and to a synthetic non-Gaussian surrogate consistent with the same data. A summary and conclusion are provided in the last section.

#### 2. The Method

##### 2.1. Cumulant-Based Independent Subspace Analysis

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

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

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

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

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

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

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

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

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

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

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

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

##### 2.2. Independent Subspace Analysis from Finite Samples

###### 2.2.1. Statistical Tests of Multivariate Non-Gaussianity

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

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

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

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

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