Spectral data often contains hundreds or thousands of wavelength channels, many of which are redundant or irrelevant for a given analytical task. Feature wavelength selection aims to identify a minimal subset of informative wavelengths that preserve predictive power while reducing dimensionality, noise, and computational cost.
Principal Component Analysis (PCA)
PCA transforms the original spectral space into a new coordinate system where axes (principal components) are ordered by the amount of variance they explain. Although PCA does not directly select individual wavelentghs, it provides a compact representation useful for subsequent modeling.
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
def simulate_spectra(n_samples=100, n_bands=200):
np.random.seed(42)
return np.random.rand(n_samples, n_bands)
def run_pca(data, n_comp=2):
transformer = PCA(n_components=n_comp)
transformed = transformer.fit_transform(data)
variance_ratios = transformer.explained_variance_ratio_
return transformed, variance_ratios
# Example usage
spectra = simulate_spectra()
proj, ratios = run_pca(spectra)
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.scatter(proj[:, 0], proj[:, 1], alpha=0.7)
plt.xlabel('PC1'); plt.ylabel('PC2')
plt.title('First Two Principal Components')
plt.subplot(1, 2, 2)
plt.bar(range(len(ratios)), ratios)
plt.xlabel('Component'); plt.ylabel('Explained Variance Ratio')
plt.title('Variance Explained by Each PC')
plt.tight_layout()
plt.show()
Synergy Interval Partial Least Squares (SiPLS)
SiPLS divides the full spectrum into contiguous intervals, evaluates each interval’s predictive ability using PLS regression with cross-validation, and selects top-performing intervals for model building.
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cross_decomposition import PLSRegression
from sklearn.model_selection import cross_val_score
def generate_synthetic_data(n_samples=100, n_wavelengths=200, noise=0.1):
np.random.seed(42)
wl = np.linspace(400, 2500, n_wavelengths)
X = np.random.randn(n_samples, n_wavelengths) + 10 * np.sin(wl)
y = 0.5 * X @ np.sin(wl) + noise * np.random.randn(n_samples)
return X, y, wl
def sipls_interval_scores(X, y, wavelengths, interval_len=10):
n_int = X.shape[1] // interval_len
errors = []
for i in range(n_int):
start, end = i * interval_len, (i + 1) * interval_len
model = PLSRegression(n_components=2)
mse = -cross_val_score(model, X[:, start:end], y, cv=5,
scoring='neg_mean_squared_error').mean()
errors.append(mse)
centers = wavelengths[interval_len//2 :: interval_len][:n_int]
return np.array(errors), centers
X, y, wl = generate_synthetic_data()
scores, centers = sipls_interval_scores(X, y, wl)
plt.figure(figsize=(12, 6))
plt.subplot(2, 1, 1)
for spec in X:
plt.plot(wl, spec, color='gray', alpha=0.3)
plt.title('Simulated Spectra')
plt.subplot(2, 1, 2)
plt.plot(centers, scores, 'o-')
plt.xlabel('Interval Center Wavelength (nm)')
plt.ylabel('Cross-Validated MSE')
plt.title('SiPLS Interval Performance')
plt.tight_layout()
plt.show()
Successive Projections Algorithm (SPA)
SPA iteratively selects variibles that are maximally orthogonal to those already chosen, minimizing collinearity while preserving information relevant to the response.
import numpy as np
import matplotlib.pyplot as plt
def spa_select(X, k):
m, n = X.shape
selected = []
P = np.eye(m)
for _ in range(k):
proj = X.T @ P
norms = np.sum(proj**2, axis=0)
idx = np.argmax(norms)
selected.append(idx)
xi = X[:, [idx]]
denom = xi.T @ P @ xi
if denom != 0:
P -= (P @ xi @ xi.T @ P) / denom
return selected
np.random.seed(42)
X_sim = np.random.rand(200, 100)
wavelengths = np.arange(100)
chosen = spa_select(X_sim, 10)
avg_spec = X_sim.mean(axis=0)
plt.figure(figsize=(10, 5))
plt.plot(wavelengths, avg_spec, label='Mean Spectrum')
plt.scatter(wavelengths[chosen], avg_spec[chosen], c='red', label='Selected by SPA')
plt.xlabel('Wavelength Index'); plt.ylabel('Intensity')
plt.legend(); plt.title('SPA-Selected Features')
plt.show()
Competitive Adaptive Reweighted Sampling (CARS)
While a full CARS implementation involves exponential decay and Monte Carlo sampling, a practical approximation uses univariate statistical tests to rank and select features.
import numpy as np
import matplotlib.pyplot as plt
from sklearn.feature_selection import SelectKBest, f_regression
np.random.seed(42)
X_raw = np.random.rand(100, 200)
y_target = np.random.rand(100)
selector = SelectKBest(score_func=f_regression, k=10)
X_reduced = selector.fit_transform(X_raw, y_target)
mask = selector.get_support()
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.imshow(X_raw.T, aspect='auto', cmap='plasma')
plt.title('Original Spectra')
plt.colorbar()
plt.subplot(1, 2, 2)
plt.imshow(X_reduced.T, aspect='auto', cmap='plasma')
plt.title('Top 10 Features (F-test)')
plt.colorbar()
plt.tight_layout()
plt.show()
Random Frog Algorithm
Random Frog approximates feature importance by repeatedly training models on random subsets and acccumulating feature weights. Here, a simplified version uses Random Forest’s built-in importance.
import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestClassifier
np.random.seed(42)
X_demo = np.random.rand(100, 200)
y_demo = np.random.randint(0, 2, 100)
def random_frog_simple(X, y, k=20, iterations=100):
importance = np.zeros(X.shape[1])
for _ in range(iterations):
idxs = np.random.choice(X.shape[1], k, replace=False)
model = RandomForestClassifier(n_estimators=10, random_state=0)
model.fit(X[:, idxs], y)
importance[idxs] += model.feature_importances_
final = np.argsort(importance)[-k:]
return np.sort(final)
selected_idx = random_frog_simple(X_demo, y_demo, k=20)
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.imshow(X_demo.T, aspect='auto', cmap='viridis')
plt.title('Original Data')
plt.subplot(1, 2, 2)
plt.imshow(X_demo[:, selected_idx].T, aspect='auto', cmap='viridis')
plt.title('Random Frog Selected Features')
plt.tight_layout()
plt.show()