2.12 Dimensionality Reduction

Geoscience datasets often carry many correlated features: oxide concentrations in a rock analysis, thousands of grid cells in a climate field, dozens of waveform attributes. Reducing the dimensionality before modeling helps because:

1. The cost of most algorithms grows with the number of input dimensions.
2. Redundant features add computation without adding information.
3. Simpler models are more robust on small datasets.
4. Fewer features make the data easier to understand.
5. Visualization is easier in two or three dimensions.

Dimensionality reduction techniques fall into two categories: feature selection and feature extraction.

1. Feature Selection¶

Feature selection keeps a subset of the original dimensions. A forward selection approach starts with the single variable that reduces the error the most and adds variables one by one. A backward selection starts with all variables and removes them one by one.

A quick first step is to look at the correlation matrix: strongly correlated features carry redundant information, and one of them can often be dropped.

We use a synthetic geochemical table from the course package mlgeo_synth. Each row is a whole-rock analysis: seven major-element oxides in wt%, density, magnetic susceptibility, and a lithology label (basalt, andesite, or granite).

🖥️ Lecture slides — Session 10 (Wed Oct 21)

SiO2 is strongly anticorrelated with MgO, FeO, and CaO, and density tracks the mafic oxides. Two effects drive this structure. First, igneous differentiation: as a melt evolves, SiO2 and K2O rise while MgO, FeO, and CaO fall. Second, closure: the oxides sum to roughly 100 wt%, so if one goes up, the others must come down. Which features would you drop based on this matrix?

2. Feature Extraction¶

Feature extraction builds a new, smaller set of dimensions as combinations of the original ones. Methods can be unsupervised (Principal Component Analysis, Independent Component Analysis) or supervised (Linear Discriminant Analysis).

3. Principal Component Analysis¶

PCA is an unsupervised method that maps the data to a lower dimensional space with minimum loss of variance.

Let $\mathbf{Y} = \mathbf{y}_1,\cdots,\mathbf{y}_n$ be the data, measured $n$ times over multiple fields of measurements (the length of $\mathbf{y}$). Each column of $\mathbf{Y}$ represents a unique observation. Each row of $\mathbf{Y}$ represents a single parameter.

To perform PCA:

1. Center the data by subtracting the mean of each row of $\mathbf{Y}$ (and usually scale each row to unit variance).
2. Calculate the covariance matrix of the centered data, $\mathbf{C} = \frac{1}{n-1} \mathbf{Y}^{\ast}\mathbf{Y}$. The covariance matrix is symmetric positive semi-definite, so it can be diagonalized.
3. Calculate the Singular Value Decomposition (SVD):

$\mathbf{X} = \mathbf{U} \Sigma \mathbf{V}^T,$

where the columns of $\mathbf{V}$ are the eigenvectors, or principal components. The first principal component points in the direction of highest variance.

3.1 The geometry of PCA: a rotated Gaussian cloud¶

To build intuition, we start with a two-dimensional point cloud: 10,000 observations drawn from a stretched, rotated Gaussian.

Step 1: subtract the mean¶

Step 2: SVD of the covariance matrix¶

The eigenvalues are close to $\sigma^2 = [4, 0.25]$, the squared lengths of the axes we used to build the cloud.

Step 3: explore the outcome¶

The projection rotates the cloud so that the direction of largest variance lies along the horizontal axis. PCA found the rotation we used to generate the data.

3.2 PCA on a Geochemical Table¶

Now we apply PCA to the geochemical table from Section 1. Features must be standardized first: oxides span tens of wt% while magnetic susceptibility is of order 10^-3 SI, and without scaling the large-magnitude features would dominate the covariance.

One component captures most of the variance, and two capture nearly all of it. The data live on a much lower dimensional surface than the nine measured features suggest.

The loadings tell us what the components mean. Each principal component is a weighted combination of the original features; the weights are called loadings.

On PC1, SiO2, K2O, and Na2O load with one sign while MgO, FeO, CaO, density, and magnetic susceptibility load with the other. This is exactly the correlation structure we saw in Section 1: oxide closure plus igneous differentiation. PC1 acts as a differentiation index. A sample’s PC1 score tells you where it sits on the basalt-to-granite spectrum, in a single number.

Note that the sign of a component is arbitrary: the SVD can return either orientation, so only the relative signs of the loadings matter.

PCA never saw the labels, yet the three lithologies separate along PC1 because composition and lithology are driven by the same underlying process. This is a common and useful outcome: an unsupervised method recovers a physically meaningful axis.

A practical note: full SVD is expensive for large matrices. Scikit-learn switches to a randomized PCA solver automatically when the data are larger than 500 x 500 and the number of requested components is less than 80% of the smaller dimension.

3.3 EOF Analysis of a Climate Field¶

Applied to spatio-temporal data, PCA yields two linked objects:

• Empirical Orthogonal Functions (EOFs): the spatial eigenvectors of the data covariance. Each EOF is a map that explains a portion of the total variance. In climate science, EOFs identify dominant patterns such as circulation modes or temperature anomaly structures.
• Principal Components (PCs): the time series that says how strongly each EOF is expressed at each time step.

Together, EOFs and PCs describe the spatial-temporal variability of the dataset.

We use mlgeo_synth.climate_field, which generates 30 years of monthly temperature anomalies on a global grid. The generator plants known structures — a seasonal mode, a zonal (land/ocean-like) mode, and a warming trend — and returns them in a truth dictionary, so we can check whether EOF analysis recovers them.

Area weighting. The grid is equal-angle: cells are spaced evenly in latitude and longitude. But the physical area of a cell shrinks toward the poles as $\cos(\phi)$. Without correction, the covariance matrix over-represents high latitudes — many grid cells, little actual area. The standard fix is to multiply each grid point by $\sqrt{\cos(\phi)}$ before the SVD, so that each cell’s contribution to the variance (which is quadratic in the data) is proportional to its area.

Did we recover the planted structure? The truth dictionary contains the seasonal and zonal patterns the generator used. We compare them to the recovered EOFs with a spatial pattern correlation. The sign of an EOF is arbitrary (an EOF and its negative describe the same mode, with the PC flipped to match), so we look at the magnitude of the correlation.

EOF1 matches the planted seasonal pattern and EOF2 matches the planted zonal pattern, with correlations at +/-1. A correlation of -1 is as good as +1 here: it is the same mode with the map and its PC both flipped. Look also at the PC time series: the seasonal PC oscillates with a 12-month period, the zonal PC varies without a trend, and PC3 — whose map is concentrated at high northern latitudes — drifts steadily in one direction. That is the planted warming trend of 0.25 deg C per decade, amplified toward the Arctic (whether the drift appears upward or downward again depends on the arbitrary sign of the EOF).

Limitations of PCA on spatio-temporal data. EOFs are constrained to be orthogonal, but physical modes of variability are not, so a single EOF can mix several processes and split others. PCA is linear, so nonlinear dynamics spread across many components. Large-scale trends can dominate the leading modes and hide local signals. And the results are sensitive to preprocessing choices: whether to remove the seasonal cycle, how to scale variables, and how to weight the grid. Treat EOFs as a description of variance, not automatically as physical modes.

4. Independent Component Analysis¶

Independent Component Analysis (ICA) separates a multivariate signal into additive, statistically independent, non-Gaussian components. It is a form of blind source separation.

Differences from PCA:

• PCA finds orthogonal axes that maximize variance, using second-order statistics (covariance). Its components are uncorrelated but not necessarily independent.
• ICA finds statistically independent components, not necessarily orthogonal, by exploiting non-Gaussianity. It requires the sources to be non-Gaussian.

In the geosciences, ICA is used for blind source separation when several unknown processes are mixed in the measurements — for example, separating earthquake, hydrologic, and seasonal contributions in geodetic time series.

The classic demonstration: three known source signals are mixed into three “receivers”, and FastICA unmixes them.

ICA recovers the three sources (up to order, sign, and scale). PCA does not: its orthogonal, maximum-variance components remain mixtures of the sources.

5. t-SNE for Visualization¶

PCA is linear. t-distributed Stochastic Neighbor Embedding (t-SNE) is a nonlinear method built for visualization: it places points in 2D so that neighbors in the high-dimensional space stay neighbors in the plane. It preserves local structure well, but distances between clusters in a t-SNE plot are not meaningful, and it is too slow for large datasets — so we subsample.

The perplexity parameter sets roughly how many neighbors each point considers. Small values fragment the data into many small clumps; large values blur local detail. Always try a few values.

Both embeddings separate the lithologies, but the geometry changes with perplexity — a reminder that t-SNE plots are qualitative. UMAP is a popular, faster alternative with similar goals; it is not installed in the course environment, but you can add it with the umap-learn package if you want to compare.

6. Other Techniques¶

1. Random projections
2. Multidimensional scaling
3. Isomap
4. Linear discriminant analysis (supervised)