Lesson 6.1.2: Neural networks for spectral classification

flowchart TD
    A[Observed Spectrum] --> B[Pre‑processing (normalisation, noise filtering)]
    B --> C[Feature Extraction (e.g. PCA, line indices)]
    C --> D[Input Layer → Hidden Layer(s)]
    D --> E[Output Layer (class probabilities)]
    E --> F[Spectral Class Assignment]
    F --> G[Evaluation (confusion matrix, accuracy)]

Concept: Foundations of Neural Networks for Spectral Classification

Stellar spectra encode a wealth of physical information: effective temperature, surface gravity, chemical composition, radial velocity and more. Extracting this information efficiently requires a framework that can learn complex, non‑linear mappings from high‑dimensional data to discrete or continuous labels. Artificial neural networks (ANNs) provide exactly this capability. In this section we unpack the theoretical underpinnings of neural networks, explain how they are adapted to spectral data, and illustrate the workflow with a minimal Python example.

1. The perceptron: the building block

A perceptron is a simple computational unit that receives a vector of inputs \( \mathbf{x} = (x_1,\dots,x_n) \), multiplies each by a weight \( w_i \), adds a bias \( b \), and applies an activation function \( \sigma \) to produce an output:

In spectral classification, each input might be the flux at a particular wavelength pixel. The weight vector encodes how strongly that pixel influences the network’s decision. The bias shifts the activation threshold, allowing the perceptron to be more or less sensitive to the combined input.

2. From a single neuron to a multi‑layer network

A single perceptron can only model linear separable problems. To capture the richly structured patterns in spectra—absorption line blends, continuum curvature, instrumental artefacts—ANNs stack many perceptrons into layers. A typical feed‑forward architecture for spectral classification consists of:

• Input layer: one node per pixel (or per feature extracted from the spectrum).
• Hidden layers: one or more layers where neurons compute weighted sums and non‑linear activations. The number of hidden units and layers is a hyper‑parameter tuned to balance expressiveness and over‑fitting.
• Output layer: for classification, a softmax layer producing a probability distribution over spectral classes; for regression, a linear or ReLU layer predicting continuous parameters.

Training adjusts the weights and biases to minimise a loss function. For classification, categorical cross‑entropy is common; for regression, mean squared error (MSE) is typical. Gradient descent, often with a variant such as Adam, propagates error gradients back through the network (back‑propagation), updating parameters iteratively.

3. Handling spectral data: preprocessing and feature engineering

Raw spectra can be irregularly sampled, noisy, and affected by telluric absorption or instrumental response. Before feeding them into a neural network, several preprocessing steps are standard:

1. Flux normalisation: divide by the continuum to remove large‑scale variations.
2. Resampling: interpolate onto a common wavelength grid to ensure each input vector has the same dimensionality.
3. Noise filtering: apply a Savitzky–Golay filter or wavelet denoising to suppress high‑frequency noise that could mislead the network.
4. Dimensionality reduction: optional use of Principal Component Analysis (PCA) or autoencoders to compress the data while preserving variance critical for classification.

These steps are crucial because neural networks are highly sensitive to the scale and distribution of inputs. Standardising each feature to zero mean and unit variance (z‑score normalisation) is a common practice that stabilises training.

4. Why neural networks excel at spectral classification

Stellar spectra are characterised by subtle, high‑dimensional patterns that are difficult to capture with simple template matching. Neural networks offer several advantages:

• Learning non‑linear feature interactions: The hierarchical layers allow the network to learn composite features, such as the relative depth of two absorption lines, that are not linearly separable.
• Robustness to noise: With sufficient training data, a network can learn to disregard random noise while focusing on the underlying spectral signatures.
• Scalability: Once trained, inference is computationally cheap, enabling real‑time classification of large survey datasets.
• Transferability: Models trained on one instrument’s spectra can be fine‑tuned for another, leveraging shared physical features across datasets.

5. A minimal implementation example

Below is a short Python snippet using the popular torch library to define a simple feed‑forward neural network for binary spectral classification (e.g. hot vs. cool stars). The code demonstrates the essential components: model definition, loss calculation, and a single training step.

import torch
import torch.nn as nn
import torch.optim as optim

1. Define a simple network
class SpectralNet(nn.Module):
    def __init__(self, input_dim, hidden_dim):
        super(SpectralNet, self).__init__()
        self.fc1 = nn.Linear(input_dim, hidden_dim)   # Input -> Hidden
        self.relu = nn.ReLU()                         # Non‑linear activation
        self.fc2 = nn.Linear(hidden_dim, 1)           # Hidden -> Output
        self.sigmoid = nn.Sigmoid()                   # Binary output

    def forward(self, x):
        out = self.fc1(x)
        out = self.relu(out)
        out = self.fc2(out)
        out = self.sigmoid(out)
        return out

2. Hyper‑parameters
input_dim = 1024          # e.g. 1024 wavelength pixels
hidden_dim = 128
learning_rate = 0.001

3. Instantiate model, loss, optimiser
model = SpectralNet(input_dim, hidden_dim)
criterion = nn.BCELoss()               # Binary cross‑entropy
optimizer = optim.Adam(model.parameters(), lr=learning_rate)

4. Dummy training data (batch of 32 spectra)
Assume X is a tensor of shape (32, 1024) and y is (32, 1)
X = torch.randn(32, input_dim)
y = torch.randint(0, 2, (32, 1)).float()

5. Forward pass
outputs = model(X)                      # Shape: (32, 1)

6. Compute loss
loss = criterion(outputs, y)

7. Backward pass and optimisation
optimizer.zero_grad()                   # Clear gradients from previous step
loss.backward()                         # Compute gradients
optimizer.step()                        # Update parameters

print(f'Loss: {loss.item():.4f}')

Line‑by‑line explanation:

• import torch – imports the PyTorch library.
• Class SpectralNet inherits from nn.Module and defines the network architecture: a linear layer to the hidden layer, a ReLU activation, another linear layer to the output, and a sigmoid to squash the output to [0,1].
• In __init__, nn.Linear creates a weight matrix and bias vector for each layer.
• The forward method specifies how data flows through the network.
• Hyper‑parameters set the dimensionality of input and hidden layers, plus the learning rate.
• We instantiate the model, the binary cross‑entropy loss (appropriate for two‑class problems), and the Adam optimiser.
• For illustration, random tensors X and binary targets y simulate a mini‑batch of spectra.
• The forward pass returns predictions outputs.
• The loss is computed by comparing predictions to true labels.
• Gradients are cleared, back‑propagated, and the optimiser updates the model parameters.
• The loss value is printed, giving a quick diagnostic of training progress.

6. Common pitfalls and best practices

While the example above is deliberately simple, real‑world spectral classification pipelines confront several challenges:

• Class imbalance: Some spectral types may be under‑represented. Techniques such as class weighting, focal loss, or oversampling can mitigate bias.
• Over‑fitting: Deep networks with many parameters can memorise noise. Regularisation methods—dropout, L2 weight decay, early stopping—help generalise.
• Spectral resolution mismatch: Training data from one instrument may differ in resolution from test data. Convolutional layers can learn to be resolution‑invariant, or data can be resampled to a common resolution.
• Interpretability: Unlike template fitting, neural networks are often seen as black boxes. Techniques such as layer‑wise relevance propagation (LRP) or saliency maps can highlight which spectral regions drive predictions.
• Hardware constraints: Large input dimensions (high‑resolution spectra) can strain memory. Batch processing and mixed‑precision training (FP16) can alleviate this.

7. Outlook: from classification to parameter inference

Once a network reliably assigns a spectral class, the same architecture can be repurposed for regression tasks, such as estimating effective temperature or metallicity. This typically involves replacing the sigmoid output with a linear activation and using MSE loss. Transfer learning—initialising a regression network with the weights of a pre‑trained classifier—often speeds up convergence and improves performance, especially when training data are scarce.

8. Quick example: Convolutional 1‑D classifier

For spectra, local wavelength patterns can be captured with a 1‑D convolutional layer before the dense blocks. The following minimal example demonstrates how to replace the first dense layer with a convolutional feature extractor.

model = models.Sequential(name='spectral_cnn_classifier')
Convolutional feature extractor
model.add(layers.Conv1D(32, kernel_size=5, activation='relu',
                       input_shape=(X_train.shape[1], 1),
                       kernel_initializer='he_normal'))
model.add(layers.MaxPooling1D(pool_size=2))
model.add(layers.Dropout(0.3))

Flatten and dense layers
model.add(layers.Flatten())
model.add(layers.Dense(128, activation='relu',
                       kernel_initializer='he_normal'))
model.add(layers.Dropout(0.3))
model.add(layers.Dense(num_classes, activation='softmax',
                       kernel_initializer='glorot_uniform'))

Compile and train as before
optimizer = optimizers.Adam(learning_rate=0.001)
model.compile(optimizer=optimizer,
              loss='categorical_crossentropy',
              metrics=['accuracy'])
history = model.fit(X_train[..., np.newaxis], y_train,
                    epochs=200,
                    batch_size=64,
                    validation_data=(X_val[..., np.newaxis], y_val),
                    callbacks=[early_stop, checkpoint],
                    verbose=2)

Note the reshaping of inputs with [..., np.newaxis] to provide a channel dimension required by Conv1D.

Worked Example: Classifying Real‑World Spectra

In this worked example we walk through every stage of a typical neural‑network workflow for classifying astronomical spectra. The example uses the publicly available Sloan Digital Sky Survey (SDSS) dataset, specifically the spectroscopic catalogue of galaxies and quasars, to illustrate how raw data are transformed into a machine‑learning ready format, how a model is trained and evaluated, and finally how the trained model can be deployed for production inference. The steps are organised into three logical units: data preparation, model development, and deployment. Each unit is described in detail, with code snippets that can be run in a standard Python 3 environment with the usual scientific stack (NumPy, Pandas, Scikit‑Learn, TensorFlow/Keras).

1. Data Acquisition and Pre‑processing

The SDSS spectra are distributed as FITS files, each containing a flux vector, an error vector, and a wavelength array. For this example we have pre‑selected 10 000 spectra: 7 000 galaxies, 2 000 quasars, and 1 000 stars. These files are stored in a directory structure that mirrors the SDSS plate‑IFU naming convention. The first task is to read the FITS files, normalise the flux, and label each spectrum according to its class.

Because the raw flux values span several orders of magnitude, we apply a log‑transform after adding a small offset to avoid taking the logarithm of zero:

import os
import numpy as np
import pandas as pd
from astropy.io import fits

def read_spectrum(file_path):
    with fits.open(file_path) as hdul:
        flux = hdul[1].data['flux']      # 1‑D array of flux values
        error = hdul[1].data['error']    # corresponding errors
        wave = hdul[1].data['wavelength']  # wavelength grid
    return wave, flux, error

def log_normalise(flux, eps=1e-10):
    return np.log10(flux + eps)

def load_dataset(root_dir):
    spectra = []
    labels = []
    for cls in ['galaxy', 'quasar', 'star']:
        cls_dir = os.path.join(root_dir, cls)
        for fname in os.listdir(cls_dir):
            if fname.endswith('.fits'):
                wave, flux, error = read_spectrum(os.path.join(cls_dir, fname))
                spectra.append(log_normalise(flux))
                labels.append(cls)
    return np.array(spectra), np.array(labels)

Line‑by‑line explanation:

• import os – standard library for file system access.
• import numpy as np, import pandas as pd – numerical and tabular tools.
• from astropy.io import fits – FITS file reader.
• def read_spectrum(file_path) – opens a FITS file, pulls the flux, error, and wavelength arrays from the second HDU (HDU index 1).
• def log_normalise(flux, eps=1e-10) – applies a base‑10 logarithm after adding a tiny epsilon to avoid log(0).
• def load_dataset(root_dir) – iterates over the three class subdirectories, reads each spectrum, normalises it, and stores the label.
• np.array(spectra) – converts the list of arrays into a 2‑D NumPy array where each row is a spectrum.
• Labels are stored as a NumPy array of strings.

After loading, we perform a train–validation–test split using Scikit‑Learn’s train_test_split, ensuring stratified sampling so that each class is represented proportionally in all subsets.

from sklearn.model_selection import train_test_split

X, y = load_dataset('data/sdss_spectra')
X_train, X_temp, y_train, y_temp = train_test_split(
    X, y, test_size=0.3, stratify=y, random_state=42)
X_val, X_test, y_val, y_test = train_test_split(
    X_temp, y_temp, test_size=0.5, stratify=y_temp, random_state=42)

We then apply a standardisation pipeline that centres each spectrum to zero mean and unit variance. This is essential because the neural network is sensitive to the scale of its inputs.

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_val_scaled = scaler.transform(X_val)
X_test_scaled = scaler.transform(X_test)

2. Model Development and Evaluation

For spectral classification a one‑dimensional convolutional neural network (1‑D CNN) is well suited, as it can learn local spectral features such as emission or absorption lines. The following architecture consists of two convolutional blocks followed by a dense classifier.

import tensorflow as tf
from tensorflow.keras import layers, models, optimizers, callbacks

def build_cnn(input_shape, num_classes):
    model = models.Sequential()
    # First convolutional block
    model.add(layers.Conv1D(32, kernel_size=5, activation='relu',
                            input_shape=input_shape))
    model.add(layers.BatchNormalization())
    model.add(layers.MaxPooling1D(pool_size=2))
    # Second convolutional block
    model.add(layers.Conv1D(64, kernel_size=5, activation='relu'))
    model.add(layers.BatchNormalization())
    model.add(layers.MaxPooling1D(pool_size=2))
    # Flatten and dense head
    model.add(layers.Flatten())
    model.add(layers.Dense(128, activation='relu'))
    model.add(layers.Dropout(0.5))
    model.add(layers.Dense(num_classes, activation='softmax'))
    return model

Line‑by‑line explanation:

• model = models.Sequential() – creates a linear stack of layers.
• Conv1D(32, kernel_size=5, activation='relu', input_shape=input_shape) – 32 filters of width 5, ReLU non‑linearity, expects input shape (n_steps, 1).
• BatchNormalization() – stabilises training by normalising activations.
• MaxPooling1D(pool_size=2) – reduces the temporal dimension by a factor of two.
• Second block doubles the filter count to 64, enabling more complex feature extraction.
• Flatten() – converts the feature map into a 1‑D vector.
• Dense(128, activation='relu') – fully connected layer with 128 units.
• Dropout(0.5) – randomly zeroes 50 % of the units during training to prevent over‑fitting.
• Dense(num_classes, activation='softmax') – final classifier; softmax ensures the outputs sum to one.

We compile the model with categorical cross‑entropy loss and the Adam optimiser, and we include a learning‑rate scheduler and early stopping callbacks to avoid unnecessary epochs.

num_classes = len(np.unique(y_train))
model = build_cnn((X_train_scaled.shape[1], 1), num_classes)

model.compile(optimizer=optimizers.Adam(learning_rate=1e-3),
              loss='sparse_categorical_crossentropy',
              metrics=['accuracy'])

lr_scheduler = callbacks.ReduceLROnPlateau(
    monitor='val_loss', factor=0.5, patience=3, verbose=1)

early_stop = callbacks.EarlyStopping(
    monitor='val_loss', patience=6, restore_best_weights=True, verbose=1)

history = model.fit(
    X_train_scaled[..., np.newaxis],  # add channel dimension
    y_train,
    validation_data=(X_val_scaled[..., np.newaxis], y_val),
    epochs=50,
    batch_size=64,
    callbacks=[lr_scheduler, early_stop],
    verbose=2)

After training, we evaluate the model on the held‑out test set. The confusion matrix gives a clear picture of class‑wise performance.

from sklearn.metrics import classification_report, confusion_matrix
import matplotlib.pyplot as plt
import seaborn as sns

y_pred = model.predict(X_test_scaled[..., np.newaxis])
y_pred_classes = np.argmax(y_pred, axis=1)

print(classification_report(y_test, y_pred_classes,
                            target_names=['galaxy', 'quasar', 'star']))

cm = confusion_matrix(y_test, y_pred_classes)
plt.figure(figsize=(6,5))
sns.heatmap(cm, annot=True, fmt='d',
            xticklabels=['galaxy', 'quasar', 'star'],
            yticklabels=['galaxy', 'quasar', 'star'],
            cmap='Blues')
plt.xlabel('Predicted class')
plt.ylabel('True class')
plt.title('Confusion matrix')
plt.show()

The report typically shows an accuracy above 95 % for galaxies and quasars, with a slightly lower rate for stars due to their spectral similarity to some late‑type galaxies. The confusion matrix visualises the misclassifications, revealing that most stellar spectra are confused with dwarf galaxies.

2.1 Alternative Architecture: Fully‑Connected MLP

As an alternative to convolutional layers, a simple multi‑layer perceptron (MLP) can be used when the dataset is small or when computational resources are limited. The MLP treats each wavelength bin as an independent feature, learning global relationships via dense layers.

def build_mlp(input_shape, num_classes):
    model = models.Sequential()
    model.add(layers.Flatten(input_shape=input_shape))
    model.add(layers.Dense(256, activation='relu'))
    model.add(layers.BatchNormalization())
    model.add(layers.Dropout(0.3))
    model.add(layers.Dense(128, activation='relu'))
    model.add(layers.BatchNormalization())
    model.add(layers.Dropout(0.3))
    model.add(layers.Dense(num_classes, activation='softmax'))
    return model

mlp_model = build_mlp((X_train_scaled.shape[1], 1), num_classes)
mlp_model.compile(optimizer=optimizers.Adam(1e-3),
                  loss='sparse_categorical_crossentropy',
                  metrics=['accuracy'])
mlp_history = mlp_model.fit(
    X_train_scaled[..., np.newaxis], y_train,
    validation_data=(X_val_scaled[..., np.newaxis], y_val),
    epochs=50, batch_size=64,
    callbacks=[lr_scheduler, early_stop], verbose=2)

Although the MLP may be slightly slower to train, it can achieve comparable performance on this dataset, especially when the spectral resolution is modest.

3. Model Deployment

Once a satisfactory model has been achieved, the next step is to expose it as an inference service. In a research environment we might use a simple Flask API that accepts a raw spectrum, applies the same pre‑processing, and returns the predicted class. For production settings, a TensorFlow SavedModel format can be exported and served via TensorFlow Serving or a cloud inference endpoint.

import joblib
import tensorflow as tf

Persist the scaler
joblib.dump(scaler, 'models/scaler.pkl')

Save the model in TensorFlow SavedModel format
model.save('models/spectral_classifier')

The Flask wrapper demonstrates how the pipeline would be applied to new data:

from flask import Flask, request, jsonify
import numpy as np
import joblib
import tensorflow as tf

app = Flask(__name__)
scaler = joblib.load('models/scaler.pkl')
model = tf.keras.models.load_model('models/spectral_classifier')

@app.route('/predict', methods=['POST'])
def predict():
    # Expect JSON with a key 'flux' containing a list of values
    data = request.get_json()
    flux = np.array(data['flux'])
    flux_log = np.log10(flux + 1e-10)
    flux_scaled = scaler.transform([flux_log])
    pred = model.predict(flux_scaled[..., np.newaxis])
    class_idx = np.argmax(pred, axis=1)[0]
    classes = ['galaxy', 'quasar', 'star']
    return jsonify({'prediction': classes[class_idx]})

if __name__ == '__main__':
    app.run(host='0.0.0.0', port=5000)

4. Performance Monitoring and Re‑training

Deploying a model is not a one‑off event. Spectral quality can change over time due to instrument degradation, and new classes (e.g. rare exotic objects) may appear in the data stream. Hence, continuous evaluation is essential. A simple strategy is to log incoming spectra and their predictions, then periodically re‑run the training pipeline on a curated subset that includes any misclassified examples. This incremental learning approach keeps the model robust without retraining from scratch.

In addition, metrics such as the Receiver Operating Characteristic (ROC) curve for each class can be plotted over successive deployments to detect drift. If the area under the curve (AUC) falls below a threshold, a full retraining cycle is triggered.

By following this end‑to‑end workflow—starting from raw FITS files, through rigorous pre‑processing, a carefully designed 1‑D CNN, thorough evaluation