The technique of Empirical Orthogonal Function analysis, usually just referred to as EOF analysis, is very common in geophysical sciences. The general aim of EOF analysis is to simplify a spatial-temporal data set by transforming it to spatial patterns of variability and temporal projections of these patterns. The spatial patterns are the EOFs, and can be thought of as basis functions in terms of variance. The associated temporal projections are the pricipal components (PCs) and are the temporal coefficients of the EOF patterns.

EOF analysis is really just re-expressing the original data set in terms of a variance basis. The original data set can be completely reconstructed using the EOFs and PCs. iHowever, in practice it is often only a subset of the EOFs that are of interest. Individual EOFs can sometimes have a physical interpretation assigned to them, or sometimes one wishes to produce a data set that is truncated in terms of variance by reconstructing the original data set using only a limited number of EOFs.

A method based on singular value decomposition (SVD) is used in `eof2` [1]. This avoids having to compute the covariance matrix directly and is therefore optimal for data sets with a large spatial dimension.

The input to EOF analysis is a spatial-temporal field. This is represented in Python by an array (or `cdms2` variable) of two or more dimensions. When one inputs an array or variable into an `eof2` solver it is re-shaped and stored internally as a two-dimensional array where the first dimension is time and the second dimension is space. It is a formal requirement of EOF analysis that this array have columns with zero-mean. The `eof2` solvers will automatically attempt to subtract the mean from each column unless the keyword argument *center* is set to *False*. This may be useful when anomalies have already been computed.

Any missing values in the two-dimensional data array are identified and removed [2]. The SVD of the anomaly matrix is then computed. The SVD computed is the truncated form, where only singular vectors (EOFs/PCs) that correspond to non-zero singular values are returned. Since the singular value is related to the fraction of variance represented by an EOF, neglecting those with singular values of zero retains a full solution. This method further reduces the computational cost of the analysis [3]. The EOFs are the right singular vectors and the standardized PCs are the left singular vectors while the squared singular values are the variances associated with each EOF.

Consider a data set that consists of observations of a single geophysical variable at multiple positions in space \(x_1, x_2, \ldots, x_M\) and at multiple times \(t_1, t_2, \ldots, t_N,\). These observations are arranged in a matrix \(\mathbf{F}\) with dimension \(N \times M\) such that each row of \(\mathbf{F}\) is a map of observations at all points in space at a particular time, and each column is a time-series of observations at a particular point at all times. The time-mean is then removed from of the \(M\) time series to form the anomaly matrix \(\mathbf{A}\) whose columns have zero-mean:

\[\begin{split}\mathbf{A} = \begin{pmatrix}
a_{1,1} & a_{1,2} & \cdots & a_{1,M} \\
a_{2,1} & a_{2,2} & \cdots & a_{2,M} \\
\vdots & \vdots & \ddots & \vdots \\
a_{N,1} & a_{N,2} & \cdots & a_{N,M}
\end{pmatrix}\end{split}\]

Typically one would then compute the covariance matrix \(\mathbf{R} = \mathbf{A}^\intercal \mathbf{A}\) and solve the eigenvalue problem:

(1)\[\mathbf{R C} = \mathbf{C} \Lambda\]

where the columns of \(\mathbf{C}\) are the eigenvectors (EOFs) and the eigenvalues (EOF variances) are on the leading diagonal of \(\Lambda\). The PCs \(\mathbf{P}\) can then be computed from the projection of \(\mathbf{A}\) onto the EOFs:

\[\mathbf{P} = \mathbf{A C}\]

The SVD method used by the `eof2` solvers is slightly different to this. Instead of computing the covariance matrix directly it computes the SVD of \(\mathbf{A}\):

\[\mathrm{SVD}\left(\mathbf{A}\right) = \mathbf{U} \Gamma \mathbf{V}^\intercal\]

The columns of \(\mathbf{U}\) and \(\mathbf{V}\) are the singular vectors and the singular values are on the leading diagonal of \(\Gamma\). To show relation between the SVD and (1) we express the covariance matrix in two ways, first a rearranged form of (1) [4]:

(2)\[\mathbf{R} = \mathbf{C} \Lambda \mathbf{C}^\intercal\]

and second the expression for the covariance matrix \(\mathbf{R}\) after first taking the SVD of \(\mathbf{A}\):

(3)\[\mathbf{R} = \mathbf{A}^\intercal \mathbf{A} = \left( \mathbf{U} \Gamma \mathbf{V}^\intercal \right)^\intercal \left( \mathbf{U} \Gamma \mathbf{V}^\intercal \right) = \mathbf{V} \Gamma^\intercal \mathbf{U}^\intercal \mathbf{U} \Gamma \mathbf{V}^\intercal = \mathbf{V} \Gamma^\intercal \Gamma \mathbf{V}^\intercal\]

It should be clear from (2) and (3) that \(\mathbf{C} = \mathbf{V}\) and \(\Lambda = \Gamma^\intercal \Gamma\). An extra benefit of the SVD method is that the singular vectors in \(\mathbf{U}\) are the standardized PCs. This can be shown by first forming an expression for \(\mathbf{A}\) in terms of the EOFs and the PCs, note that this is the reconstruction expression:

\[\mathbf{A} = \mathbf{P} \mathbf{C}^\intercal\]

By defining a normalized PC \(\phi_j\) as

\[\phi_j = \dfrac{\mathbf{p}_j}{\sqrt{\lambda_j}}\]

where \(\mathbf{p}_j\) is a column of \(\mathbf{P}\), and defining a diagonal matrix \(\mathbf{D}\) with \(\sqrt{\lambda}_j\) on the leading diagonal and a matrix \(\Phi\) of the ordered \(\phi_j\) as columns we can write a new expression for \(\mathbf{A}\):

\[\mathbf{A} = \Phi \mathbf{D} \mathbf{C}^\intercal\]

It should be clear that this expression is equivalent to the SVD of \(\mathbf{A}\) and therefore that the left singular vectors are the PCs scaled to unit variance.

Footnotes

[1] | Note that this is the linear algebra technique and is different to the coupled SVD analysis of two fields. The latter technique is often (and less ambiguously) referred to as maximum covariance analysis (MCA). |

[2] | An often used alternative and equally valid strategy is to set missing values to a constant value. This method yields the same solution since a constant time-series has zero-variance, but increases the computational cost of computing the EOFs. |

[3] | For an \(N\) by \(M\) anomaly matrix the rank of the corresponding covariance matrix can be at most \(\mathrm{min}\left(m,n\right)\). The number of zero eigenvalues is at least \(\left|m - n\right|\). |

[4] | This rearrangement is possible since the column eigenvectors in \(\mathbf{C}\) are mutually orthogonal and hence \(\mathbf{C} \mathbf{C}^\intercal = \mathbf{I}\) where \(\mathbf{I}\) is the identity matrix. |