LSTM for Multivariate Time Series

Forecasting a single number is hard. Forecasting when your data has several interlocking signals — sensor readings, web traffic, and sales all moving together — is harder still. A plain neural network or a standard LSTM sees each time series in isolation and ignores the fact that rising traffic usually precedes rising sales. This lesson removes that bottleneck: you will learn LSTM for multivariate time series, a technique that lets one recurrent model ingest many input channels at once and learn their joint dynamics. By the end, you will not only understand the architecture but also build a working multivariate forecaster in PyTorch, step by step.

The problem this lesson solves

Most real-world forecasting problems do not arrive as a single clean column of numbers. A factory floor sensor records temperature, vibration, and pressure simultaneously. A marketing dashboard logs ad spend, clicks, and conversions side by side. Each of those channels influences the others, and the future of one depends on the history of all.

Here is the painful scenario this lesson fixes: you have a multivariate dataset ready for machine learning, but your model keeps treating it as a set of independent univariate problems. You flatten the features, train an LSTM, and get mediocre accuracy. You tweak hyperparameters, add layers, and still the validation loss refuses to drop below a plateau. The root cause is almost never model capacity — it is input shape. A standard LSTM expects 3D input shaped as (samples, timesteps, features). If you pass (samples, features) or (samples, timesteps, 1) by accident, you are telling the model to forget everything about cross-feature dependencies.

This lesson teaches you how to structure data and build the model so that LSTM for multivariate time series works as intended. No more throwing away the relationships between your channels.

Why learn this now? If you are following the Applied AI engineering path, you have already mastered sequence basics. Multivariate input is the natural next step before you move to attention-based models and transformers.

Core concept / mental model

Think of an LSTM as a team of listeners at a meeting. Each listener (a hidden unit) pays attention to the conversation but takes notes differently. In a univariate LSTM, all listeners hear the same single speaker — the same value over time. In a multivariate LSTM, each listener hears a different speaker: one hears temperature, another hears vibration, a third hears pressure. Crucially, they sit in the same room and share notes at every step. The internal gates — forget, input, output — decide which notes to keep, which to update, and which to whisper to the next room.

That sharing is the magic. At timestep t, the hidden state does not just depend on the current value of a single feature. It depends on the entire vector of all features at time t, plus the hidden state from t−1. The LSTM learns which cross-feature patterns matter: e.g., when temperature spikes and vibration rises, pressure is likely to follow.

Key terminology to keep straight:

• Timestep: a point in the sequence (e.g., hour 12).
• Feature: one channel of input (e.g., temperature).
• Window: a sliding slice of timesteps fed to the model at once.
• Horizon: how far into the future you predict.

Here is a word-diagram of the data flow for a single window:

Input window: shape (window_size, n_features)
      |
      v
[ LSTM cell ] -- hidden state h_t, cell state c_t
      |
      v
Output: prediction for next timestep(s)

The LSTM cell consumes the whole (window_size, n_features) tensor and produces a final hidden state that summarizes what it learned. That summary feeds a dense layer to output a prediction.

How it works step by step

Getting LSTM for multivariate time series right is about preparing data and building the model in the correct shape. Follow these logical steps:

1. Gather and combine your features. All input channels must be aligned on the same time index. Missing values? Impute or drop rows. Outliers? Clip or winsorize.

2. Normalize each feature independently. LSTMs are sensitive to scale. Use MinMaxScaler or StandardScaler per channel. This prevents a large-magnitude feature like pressure from dominating the gradients.

3. Create sliding windows. For each window of W timesteps, you produce an input of shape (W, n_features) and a target of shape (1,) (or (horizon,)). This is the single most important step — if the shape is wrong, everything else collapses.

4. Split chronologically. Random shuffling destroys temporal order. Use a train/validation split that respects time (e.g., first 80% for training, last 20% for validation).

5. Feed the data as 3D tensors. The model expects (batch_size, W, n_features). Convert your arrays to torch.Tensor and add a batch dimension.

6. Build the LSTM + dense head. The LSTM processes the sequence and outputs the final hidden state. A linear layer maps that to your prediction.

7. Train with an appropriate loss. For regression, MSELoss is standard. Use an optimizer like Adam with a learning rate around 0.001.

8. Evaluate with a metric that matters. MSE is good, but RMSE and MAE are more interpretable. Plot predictions vs. actuals to see how well you track the trend.

The order matters: normalize before windowing, and split before shuffling (you won't shuffle at all).

Hands-on walkthrough

Let's build a complete, runnable example. We'll use a synthetic dataset with two related features: feature_1 (sine wave with noise) and feature_2 (feature_1 shifted plus extra noise). The target is the next value of feature_1. This mimics a real scenario where one channel leads another.

First, import libraries and create the dataset:

import numpy as np
import torch
import torch.nn as nn
from sklearn.preprocessing import MinMaxScaler

torch.manual_seed(42)
np.random.seed(42)

# Create two related time series
n = 500
t = np.arange(n)
feature_1 = np.sin(0.02 * t) + 0.3 * np.sin(0.5 * t)
feature_2 = np.roll(feature_1, 5) + 0.1 * np.random.randn(n)  # feature_1 delayed by 5 steps

# Combine into one matrix: shape (n, 2)
data = np.column_stack([feature_1, feature_2])

# Normalize each feature independently
scaler = MinMaxScaler()
data_scaled = scaler.fit_transform(data)

Next, create sliding windows. We'll use a window size of 10 timesteps and predict the next single value:

def create_sequences(data, window_size):
    X, y = [], []
    for i in range(len(data) - window_size):
        X.append(data[i:i+window_size])
        y.append(data[i+window_size, 0])  # predict feature_1's next value
    return np.array(X), np.array(y)

window_size = 10
X, y = create_sequences(data_scaled, window_size)

# Chronological split: 80% train, 20% test
split_idx = int(0.8 * len(X))
X_train, X_test = X[:split_idx], X[split_idx:]
y_train, y_test = y[:split_idx], y[split_idx:]

# Convert to tensors: (samples, window, features)
X_train_t = torch.tensor(X_train, dtype=torch.float32)
y_train_t = torch.tensor(y_train, dtype=torch.float32).unsqueeze(1)
X_test_t = torch.tensor(X_test, dtype=torch.float32)
y_test_t = torch.tensor(y_test, dtype=torch.float32).unsqueeze(1)

print(f"Train shape: {X_train_t.shape}")  # e.g., (390, 10, 2)

Now define and train the model:

class MultivariateLSTM(nn.Module):
    def __init__(self, n_features, hidden_size=32, num_layers=1):
        super().__init__()
        self.lstm = nn.LSTM(input_size=n_features, hidden_size=hidden_size, num_layers=num_layers, batch_first=True)
        self.fc = nn.Linear(hidden_size, 1)

    def forward(self, x):
        # x shape: (batch, window, n_features)
        out, _ = self.lstm(x)  # out: (batch, window, hidden_size)
        last_hidden = out[:, -1, :]  # take the last timestep's hidden state
        return self.fc(last_hidden)

model = MultivariateLSTM(n_features=2, hidden_size=32)
loss_fn = nn.MSELoss()
optimizer = torch.optim.Adam(model.parameters(), lr=0.001)

# Training loop
epochs = 100
batch_size = 32
n_batches = len(X_train_t) // batch_size

for epoch in range(epochs):
    model.train()
    total_loss = 0
    for i in range(n_batches):
        batch_X = X_train_t[i*batch_size:(i+1)*batch_size]
        batch_y = y_train_t[i*batch_size:(i+1)*batch_size]

        pred = model(batch_X)
        loss = loss_fn(pred, batch_y)

        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
        total_loss += loss.item()
    if (epoch+1) % 20 == 0:
        print(f"Epoch {epoch+1}/{epochs}, Loss: {total_loss/n_batches:.4f}")

Finally, evaluate on the test set and inspect predictions:

model.eval()
with torch.no_grad():
    predictions = model(X_test_t).numpy().flatten()
    targets = y_test_t.numpy().flatten()

# Inverse transform to original scale for interpretation
pred_inv = scaler.inverse_transform(np.column_stack([predictions, np.zeros_like(predictions)]))[:, 0]
true_inv = scaler.inverse_transform(np.column_stack([targets, np.zeros_like(targets)]))[:, 0]

print(f"First 5 predictions: {pred_inv[:5]}")
print(f"First 5 actuals:     {true_inv[:5]}")

# Quick RMSE
rmse = np.sqrt(np.mean((pred_inv - true_inv) ** 2))
print(f"Test RMSE: {rmse:.4f}")

Expected output (approximate):

Epoch 20/100, Loss: 0.0041
Epoch 40/100, Loss: 0.0028
Epoch 60/100, Loss: 0.0021
Epoch 80/100, Loss: 0.0017
Epoch 100/100, Loss: 0.0015
First 5 predictions: [-0.012, 0.104, 0.279, 0.421, 0.553]
First 5 actuals:     [-0.020, 0.098, 0.285, 0.430, 0.547]
Test RMSE: 0.0451

Your numbers will vary slightly, but the loss should steadily decrease, and predictions should track the sine wave's ups and downs.

Compare options / when to choose what

You have several ways to handle multivariate time series. Here's how they stack up:

When to choose LSTM: - Your sequences are 10–500 timesteps long. - You have enough data (thousands of samples) to train deep networks. - You need to model complex nonlinear interactions between features.

When to avoid LSTM: - You have very short sequences (< 5 timesteps) — simpler models win. - Your dataset is tiny (<100 samples) — you'll overfit. - You have strong seasonality you could model with SARIMA or Prophet.

Troubleshooting & edge cases

Even with the correct code, things can go sideways. Here are the most common issues and fixes:

Model trains but loss doesn't decrease. - Check that your data is normalized. If one feature is in thousands and another in decimals, gradients will be dominated. - Try a smaller learning rate (e.g., 0.0001). - Increase model capacity (hidden size) if you have lots of data.

Validation loss is much higher than training loss. - This is overfitting. Add regularization: dropout=0.2 in the LSTM, or weight_decay in the optimizer. - Reduce model complexity. - Increase the amount of training data.

Predictions are all the same constant. - Your targets may have low variance after normalization. Check y_train_t — if it's nearly flat, the model will predict the mean. - Your window may be too short to capture dynamics.

Shape mismatch errors. - This is the #1 error. Print X_train_t.shape and ensure it's (samples, window, n_features). If it's (samples, n_features), you forgot to create windows. - If using batch_first=False (default), the input must be (seq_len, batch, features). We set batch_first=True for convenience.

Data leakage. - Never normalize using the entire dataset's mean and variance. Fit MinMaxScaler only on the training set, then transform the test set. This simulates out-of-sample conditions.

What you learned & what's next

You now understand the core of LSTM for multivariate time series. Let's recap the key insights:

• LSTM for multivariate time series accepts input of shape (samples, timesteps, features) — each feature is a channel that the LSTM processes together.
• The sliding window approach is the bridge between raw time-stamped data and model-ready tensors.
• Normalization per feature and chronological splitting are non-negotiable for reliable forecasting.
• Your model now sees cross-feature dependencies — the hidden state at each step is a function of the entire input vector.
• You built and trained a complete PyTorch model, evaluated it with RMSE, and saw how to interpret predictions on the original scale.

You've met both learning objectives: you can explain the architecture and complete a practical hands-on exercise.

What's next? In the next lesson, you'll tackle sequence-to-sequence models — where the output isn't a single value but a whole future sequence. That builds directly on the windowing and LSTM skills you just mastered, but introduces an encoder–decoder structure. Stay with it: the same discipline you used here — careful shaping, normalization, and chronological splitting — will pay off again.