Los conjuntos de datos geocientíficos suelen cargar muchas características correlacionadas: concentraciones de óxidos en un análisis de roca, miles de celdas de malla en un campo climático, decenas de características de formas de onda. Reducir la dimensionalidad antes de modelar ayuda porque:
- El costo de la mayoría de los algoritmos crece con el número de dimensiones de entrada.
- Las características redundantes agregan cómputo sin agregar información.
- Los modelos más simples son más robustos en conjuntos de datos pequeños.
- Con menos características, los datos son más fáciles de entender.
- La visualización es más fácil en dos o tres dimensiones.
Las técnicas de reducción de dimensionalidad caen en dos categorías: la selección de características y la extracción de características.
1. Selección de características¶
La selección de características conserva un subconjunto de las dimensiones originales. Un enfoque de selección hacia adelante comienza con la única variable que más reduce el error y agrega variables una por una. Una selección hacia atrás comienza con todas las variables y las elimina una por una.
Un primer paso rápido es mirar la matriz de correlación: las características fuertemente correlacionadas cargan información redundante, y a menudo se puede descartar uno de ellos.
Usamos una tabla geoquímica sintética del paquete del curso mlgeo_synth. Cada fila es un análisis de roca total: siete óxidos de elementos mayores en % en peso, la densidad, la susceptibilidad magnética y una etiqueta de litología (basalto, andesita o granito).
# Import useful modules
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import mlgeo_synthgeochem = mlgeo_synth.geochem_table(n=5000, seed=42)
geochem.head()# The three lithology classes are imbalanced, as real map units usually are.
geochem["label"].value_counts()label
granite 2794
basalt 1713
andesite 493
Name: count, dtype: int64# Correlation matrix of the numeric features
features = geochem.drop(columns="label")
correlation_matrix = features.corr()
correlation_matrix.style.background_gradient(cmap="coolwarm", vmin=-1, vmax=1)El SiO2 está fuertemente anticorrelacionado con MgO, FeO y CaO, y la densidad sigue a los óxidos máficos. Dos efectos generan esta estructura. Primero, la diferenciación ígnea: a medida que un fundido evoluciona, SiO2 y K2O suben mientras que MgO, FeO y CaO bajan. Segundo, el cierre composicional (closure): los óxidos suman aproximadamente 100 % en peso, así que si uno sube, los demás deben bajar. ¿Qué características descartaría usted con base en esta matriz?
2. Extracción de características¶
La extracción de características construye un conjunto nuevo y más pequeño de dimensiones como combinaciones de las originales. Los métodos pueden ser no supervisados (análisis de componentes principales, análisis de componentes independientes) o supervisados (análisis discriminante lineal).
3. Análisis de componentes principales¶
El PCA es un método no supervisado que proyecta los datos a un espacio de menor dimensión con una pérdida mínima de varianza.
Sean los datos, medidos veces sobre múltiples campos de medición (la longitud de ). Cada columna de representa una observación única. Cada fila de representa un solo parámetro.
Para realizar el PCA:
- Centre los datos restando la media de cada fila de (y usualmente escale cada fila a varianza unitaria).
- Calcule la matriz de covarianza de los datos centrados, . La matriz de covarianza es simétrica y semidefinida positiva, así que puede diagonalizarse.
- Calcule la descomposición en valores singulares (SVD):
donde las columnas de son los vectores propios, o componentes principales. La primera componente principal apunta en la dirección de mayor varianza.
3.1 La geometría del PCA: una nube gaussiana rotada¶
Para construir intuición, empezamos con una nube de puntos en dos dimensiones: 10 000 observaciones extraídas de una gaussiana estirada y rotada.
# Generate the toy data
rng = np.random.default_rng(42)
xC = np.array([2, 1]) # Center of data (mean)
sig = np.array([2, 0.5]) # Principal axes
theta = np.pi / 3 # Rotate cloud by pi/3
R = np.array([[np.cos(theta), -np.sin(theta)], # Rotation matrix
[np.sin(theta), np.cos(theta)]])
nPoints = 10000
# create the cloud of points (np.matmul can also be written @)
X = R @ np.diag(sig) @ rng.standard_normal((2, nPoints)) + np.diag(xC) @ np.ones((2, nPoints))
# plot the data
fig, ax1 = plt.subplots()
ax1.plot(X[0, :], X[1, :], '.', color='k', alpha=0.125)
ax1.grid()
ax1.set_xlim((-6, 8))
ax1.set_ylim((-6, 8))
ax1.set_aspect('equal')
plt.show()
Paso 1: reste la media¶
Xavg = np.mean(X, axis=1) # Compute mean
B = X - Xavg[:, np.newaxis] # Mean-subtracted data
plt.scatter(B[0, :], B[1, :], color='k', alpha=0.125)
plt.gca().set_aspect('equal')
plt.show()
# calculate the covariance matrix
covB = (B @ B.T) / nPoints
print(f"shape of B {B.shape} and shape of covB {covB.shape}")
print(covB)shape of B (2, 10000) and shape of covB (2, 2)
[[1.21210083 1.65131023]
[1.65131023 3.08978984]]
Paso 2: SVD de la matriz de covarianza¶
U, S, VT = np.linalg.svd(covB, full_matrices=False)
print("eigenvalues (variances along each axis):", S)
print("eigenvectors (rows of VT):")
print(VT)eigenvalues (variances along each axis): [4.05048593 0.25140474]
eigenvectors (rows of VT):
[[-0.50286768 -0.8643634 ]
[-0.8643634 0.50286768]]
Los valores propios están cerca de , los cuadrados de las longitudes de los ejes que usamos para construir la nube.
Paso 3: explore el resultado¶
fig, ax2 = plt.subplots()
ax2.plot(X[0, :], X[1, :], '.', color='k', alpha=0.125) # Plot data to overlay PCA
ax2.grid()
ax2.set_xlim((-6, 8))
ax2.set_ylim((-6, 8))
ax2.set_aspect('equal')
# Plot the eigenvectors, scaled by the standard deviation along each axis
for k, color in zip(range(2), ['cyan', 'orange']):
scale = np.sqrt(S[k])
ax2.plot([Xavg[0], Xavg[0] + VT[k, 0] * scale],
[Xavg[1], Xavg[1] + VT[k, 1] * scale],
'-', color=color, linewidth=3, label=f"PC{k+1}")
ax2.legend()
plt.show()
# Project the original data onto the principal axes
projected = B.T @ VT.T
plt.scatter(projected[:, 0], projected[:, 1], c='k', alpha=0.125)
ax = plt.gca()
ax.set_axisbelow(True)
ax.grid()
ax.set_aspect('equal')
ax.set_xlabel("PC1")
ax.set_ylabel("PC2")
plt.show()
La proyección rota la nube de modo que la dirección de mayor varianza quede sobre el eje horizontal. El PCA encontró la rotación que usamos para generar los datos.
3.2 PCA sobre una tabla geoquímica¶
Ahora aplicamos el PCA a la tabla geoquímica de la sección 1. Las características deben estandarizarse primero: los óxidos abarcan decenas de % en peso mientras que la susceptibilidad magnética es del orden de 10-3 SI, y sin escalar, las características de gran magnitud dominarían la covarianza.
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
feature_names = features.columns.tolist()
scaler = StandardScaler()
geochem_scaled = scaler.fit_transform(features)
pca = PCA()
geochem_pca = pca.fit_transform(geochem_scaled)
print("Explained variance ratio:", np.round(pca.explained_variance_ratio_, 3))Explained variance ratio: [0.782 0.097 0.06 0.037 0.01 0.006 0.004 0.003 0.001]
# Scree plot: variance explained by each component
n_pc = len(pca.explained_variance_ratio_)
fig, ax = plt.subplots(figsize=(7, 4))
ax.bar(np.arange(1, n_pc + 1), pca.explained_variance_ratio_, label="per component")
ax.plot(np.arange(1, n_pc + 1), np.cumsum(pca.explained_variance_ratio_),
'o-', color='k', label="cumulative")
ax.set_xlabel("Principal component")
ax.set_ylabel("Explained variance ratio")
ax.legend()
plt.show()
Una componente captura la mayor parte de la varianza, y dos capturan casi toda. Los datos viven sobre una superficie de dimensión mucho menor de lo que sugieren las nueve características medidas.
Las cargas nos dicen qué significan las componentes. Cada componente principal es una combinación ponderada de las características originales; los pesos se llaman cargas (loadings).
fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True)
for k, ax in enumerate(axes):
loadings = pca.components_[k]
colors = ['tab:red' if v < 0 else 'tab:blue' for v in loadings]
ax.bar(feature_names, loadings, color=colors)
ax.axhline(0, color='k', linewidth=0.8)
ax.set_title(f"PC{k+1} loadings "
f"({100*pca.explained_variance_ratio_[k]:.1f}% of variance)")
ax.tick_params(axis='x', rotation=60)
axes[0].set_ylabel("Loading")
plt.tight_layout()
plt.show()
En la PC1, SiO2, K2O y Na2O cargan con un signo mientras que MgO, FeO, CaO, la densidad y la susceptibilidad magnética cargan con el otro. Esta es exactamente la estructura de correlación que vimos en la sección 1: el cierre de los óxidos más la diferenciación ígnea. La PC1 actúa como un índice de diferenciación. El puntaje de PC1 de una muestra le dice dónde se ubica en el espectro que va del basalto al granito, en un solo número.
Note que el signo de una componente es arbitrario: la SVD puede devolver cualquiera de las dos orientaciones, así que solo importan los signos relativos de las cargas.
# Scatter of the first two PCs, colored by lithology
fig, ax = plt.subplots(figsize=(8, 6))
for lith in geochem["label"].unique():
mask = (geochem["label"] == lith).to_numpy()
ax.scatter(geochem_pca[mask, 0], geochem_pca[mask, 1],
s=8, alpha=0.4, label=lith)
ax.set_xlabel("PC1 (differentiation index)")
ax.set_ylabel("PC2")
ax.legend()
ax.grid(True)
plt.show()
El PCA nunca vio las etiquetas y, sin embargo, las tres litologías se separan a lo largo de la PC1, porque la composición y la litología están gobernadas por el mismo proceso subyacente. Este es un resultado común y útil: un método no supervisado recupera un eje con significado físico.
Una nota práctica: la SVD completa es costosa para matrices grandes. Scikit-learn cambia automáticamente a un solucionador de PCA aleatorizado cuando los datos superan 500 x 500 y el número de componentes solicitadas es menor que el 80 % de la dimensión más pequeña.
3.3 Análisis EOF de un campo climático¶
Aplicado a datos espaciotemporales, el PCA produce dos objetos ligados:
- Funciones ortogonales empíricas (EOF): los vectores propios espaciales de la covarianza de los datos. Cada EOF es un mapa que explica una porción de la varianza total. En la ciencia del clima, las EOF identifican patrones dominantes como modos de circulación o estructuras de anomalías de temperatura.
- Componentes principales (PC): las series de tiempo que dicen con qué fuerza se expresa cada EOF en cada paso de tiempo.
Juntas, las EOF y las PC describen la variabilidad espacial y temporal del conjunto de datos.
Usamos mlgeo_synth.climate_field, que genera 30 años de anomalías mensuales de temperatura sobre una malla global. El generador siembra estructuras conocidas — un modo estacional, un modo zonal (tipo tierra/océano) y una tendencia de calentamiento — y las devuelve en un diccionario truth, así que podemos comprobar si el análisis EOF las recupera.
field, truth = mlgeo_synth.climate_field(
n_lat=40, n_lon=80, n_months=360, trend_c_per_decade=0.25, seed=42
)
lat = truth["lat"]
lon = truth["lon"]
n_months, n_lat, n_lon = field.shape
print("field shape (months, lat, lon):", field.shape)
print("truth keys:", list(truth.keys()))field shape (months, lat, lon): (360, 40, 80)
truth keys: ['lat', 'lon', 'seasonal_pattern', 'zonal_pattern', 'trend_c_per_decade']
# One month of the field
plt.figure(figsize=(8, 4))
plt.pcolormesh(lon, lat, field[1], cmap='coolwarm', shading='auto')
plt.title('Temperature anomaly, month 2')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.colorbar(label='deg C', fraction=0.025, pad=0.04)
plt.show()
Ponderación por área. La malla es equiangular: las celdas están espaciadas uniformemente en latitud y longitud. Pero el área física de una celda se encoge hacia los polos como . Sin corrección, la matriz de covarianza sobrerrepresenta las latitudes altas — muchas celdas de malla, poca área real. La corrección estándar es multiplicar cada punto de malla por antes de la SVD, de modo que la contribución de cada celda a la varianza (que es cuadrática en los datos) sea proporcional a su área.
# Remove the time mean at each grid point, then apply area weights
anom = field - field.mean(axis=0)
w = np.sqrt(np.cos(np.deg2rad(lat))) # shape (n_lat,)
anom_w = anom * w[None, :, None]
# Reshape to a (time x space) matrix and take the SVD
Xmat = anom_w.reshape(n_months, n_lat * n_lon)
U, S, VT = np.linalg.svd(Xmat, full_matrices=False)
variance_fraction = S**2 / np.sum(S**2)
print("Variance fraction of first 5 modes:", np.round(variance_fraction[:5], 3))Variance fraction of first 5 modes: [0.955 0.038 0.006 0. 0. ]
n_modes = 3
# Rows of VT are the EOFs of the *weighted* field; divide the weights back
# out to display physical patterns.
eofs = VT[:n_modes].reshape(n_modes, n_lat, n_lon) / w[None, :, None]
# PC time series: projection of the data on each EOF
pcs = U[:, :n_modes] * S[:n_modes]
fig, axes = plt.subplots(n_modes, 2, figsize=(12, 3 * n_modes),
gridspec_kw={'width_ratios': [1.3, 1]})
time_years = np.arange(n_months) / 12
for k in range(n_modes):
im = axes[k, 0].pcolormesh(lon, lat, eofs[k], cmap='coolwarm', shading='auto')
axes[k, 0].set_title(f"EOF{k+1} ({100*variance_fraction[k]:.1f}% of variance)")
axes[k, 0].set_ylabel("Latitude")
fig.colorbar(im, ax=axes[k, 0], fraction=0.025, pad=0.04)
axes[k, 1].plot(time_years, pcs[:, k], linewidth=0.8)
axes[k, 1].set_title(f"PC{k+1} time series")
axes[k, 1].grid(True)
axes[-1, 0].set_xlabel("Longitude")
axes[-1, 1].set_xlabel("Time (years)")
plt.tight_layout()
plt.show()
¿Recuperamos la estructura sembrada? El diccionario truth contiene los patrones estacional y zonal que usó el generador. Los comparamos con las EOF recuperadas mediante una correlación espacial de patrones. El signo de una EOF es arbitrario (una EOF y su negativo describen el mismo modo, con la PC invertida para compensar), así que miramos la magnitud de la correlación.
def pattern_corr(a, b):
"""Pearson correlation between two flattened maps."""
return np.corrcoef(a.ravel(), b.ravel())[0, 1]
planted = {"seasonal_pattern": truth["seasonal_pattern"],
"zonal_pattern": truth["zonal_pattern"]}
print(f"{'':>12s}" + "".join(f"{name:>20s}" for name in planted))
for k in range(n_modes):
row = f"{'EOF' + str(k+1):>12s}"
for name, pat in planted.items():
row += f"{pattern_corr(eofs[k], pat):>20.2f}"
print(row) seasonal_pattern zonal_pattern
EOF1 1.00 -0.00
EOF2 0.00 -1.00
EOF3 -0.12 -0.23
La EOF1 coincide con el patrón estacional sembrado y la EOF2 con el patrón zonal sembrado, con correlaciones de +/-1. Una correlación de -1 es tan buena como una de +1 aquí: es el mismo modo con el mapa y su PC ambos invertidos. Mire también las series de tiempo de las PC: la PC estacional oscila con un periodo de 12 meses, la PC zonal varía sin tendencia, y la PC3 — cuyo mapa se concentra en las latitudes altas del norte — deriva de manera sostenida en una dirección. Esa es la tendencia de calentamiento sembrada, de 0.25 grados C por década, amplificada hacia el Ártico (que la deriva aparezca hacia arriba o hacia abajo depende, otra vez, del signo arbitrario de la EOF).
Limitaciones del PCA sobre datos espaciotemporales. Las EOF están obligadas a ser ortogonales, pero los modos físicos de variabilidad no lo están, así que una sola EOF puede mezclar varios procesos y partir otros. El PCA es lineal, así que la dinámica no lineal se reparte entre muchas componentes. Las tendencias de gran escala pueden dominar los modos principales y ocultar señales locales. Y los resultados son sensibles a las decisiones de preprocesamiento: si remover o no el ciclo estacional, cómo escalar las variables y cómo ponderar la malla. Trate las EOF como una descripción de la varianza, no automáticamente como modos físicos.
4. Análisis de componentes independientes¶
El análisis de componentes independientes (ICA) separa una señal multivariada en componentes aditivas, estadísticamente independientes y no gaussianas. Es una forma de separación ciega de fuentes.
Diferencias con el PCA:
- El PCA encuentra ejes ortogonales que maximizan la varianza, usando estadísticos de segundo orden (la covarianza). Sus componentes están descorrelacionadas pero no son necesariamente independientes.
- El ICA encuentra componentes estadísticamente independientes, no necesariamente ortogonales, explotando la no gaussianidad. Requiere que las fuentes sean no gaussianas.
En las geociencias, el ICA se usa para la separación ciega de fuentes cuando varios procesos desconocidos están mezclados en las mediciones — por ejemplo, para separar las contribuciones sísmica, hidrológica y estacional en series de tiempo geodésicas.
La demostración clásica: tres señales fuente conocidas se mezclan en tres «receptores», y FastICA las desmezcla.
from scipy import signal
from sklearn.decomposition import FastICA
rng = np.random.default_rng(0)
n_samples = 2000
time = np.linspace(0, 8, n_samples)
# create 3 source signals
s1 = np.sin(2 * time) # sinusoid
s2 = np.sign(np.sin(3 * time)) # square wave
s3 = signal.sawtooth(2 * np.pi * time) # sawtooth
S_true = np.c_[s1, s2, s3]
S_true += 0.2 * rng.standard_normal(S_true.shape) # add noise
S_true /= S_true.std(axis=0) # standardize
# Mix the sources: 3 signals recorded at 3 receivers
A = np.array([[1, 1, 1], [0.5, 2, 1.0], [1.5, 1.0, 2.0]]) # mixing matrix
X_mixed = S_true @ A.T# Unmix with ICA; compare with PCA
ica = FastICA(n_components=3, random_state=0)
S_ica = ica.fit_transform(X_mixed)
pca3 = PCA(n_components=3)
S_pca = pca3.fit_transform(X_mixed)
plt.figure(figsize=(11, 8))
models = [X_mixed, S_true, S_ica, S_pca]
names = ['Observations (mixed signals)',
'True sources',
'ICA recovered signals',
'PCA recovered signals']
colors = ['red', 'steelblue', 'orange']
for ii, (model, name) in enumerate(zip(models, names), 1):
plt.subplot(4, 1, ii)
plt.title(name)
for sig_, color in zip(model.T, colors):
plt.plot(sig_, color=color)
plt.tight_layout()
plt.show()
El ICA recupera las tres fuentes (salvo el orden, el signo y la escala). El PCA no: sus componentes ortogonales de máxima varianza siguen siendo mezclas de las fuentes.
5. t-SNE para visualización¶
El PCA es lineal. El t-distributed Stochastic Neighbor Embedding (t-SNE) es un método no lineal construido para la visualización: coloca los puntos en 2D de modo que los vecinos en el espacio de alta dimensión sigan siendo vecinos en el plano. Preserva bien la estructura local, pero las distancias entre clústeres en un gráfico t-SNE no son significativas, y es demasiado lento para conjuntos de datos grandes — así que submuestreamos.
El parámetro perplexity (perplejidad) fija aproximadamente cuántos vecinos considera cada punto. Los valores pequeños fragmentan los datos en muchos grumos pequeños; los valores grandes difuminan el detalle local. Pruebe siempre varios valores.
from sklearn.manifold import TSNE
# Subsample the standardized geochemical table for speed
rng = np.random.default_rng(42)
idx = rng.choice(len(geochem), size=1500, replace=False)
X_sub = geochem_scaled[idx]
labels_sub = geochem["label"].to_numpy()[idx]
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
for ax, perp in zip(axes, [5, 50]):
emb = TSNE(n_components=2, perplexity=perp, random_state=42).fit_transform(X_sub)
for lith in np.unique(labels_sub):
mask = labels_sub == lith
ax.scatter(emb[mask, 0], emb[mask, 1], s=8, alpha=0.6, label=lith)
ax.set_title(f"t-SNE, perplexity = {perp}")
ax.set_xticks([])
ax.set_yticks([])
axes[0].legend()
plt.tight_layout()
plt.show()
Los dos embeddings (proyecciones de baja dimensión) separan las litologías, pero la geometría cambia con la perplejidad — un recordatorio de que los gráficos t-SNE son cualitativos. UMAP es una alternativa popular y más rápida con objetivos similares; no está instalada en el entorno del curso, pero puede agregarla con el paquete umap-learn si quiere comparar.
6. Otras técnicas¶
- Proyecciones aleatorias
- Escalamiento multidimensional
- Isomap
- Análisis discriminante lineal (supervisado)