N-BEATS Explained: Basis Expansion and Doubly Residual Stacking

N-BEATS is a stack of fully connected layers with two ideas bolted on: each block outputs coefficients over a basis rather than outputting the forecast directly, and each block subtracts what it explained from the input before the next one sees it. Both are visible in the shape of the tensors, and both are why it works without a single time-series-specific component.

What one block does

Boris Oreshkin, Dmitri Carpov, Nicolas Chapados and Yoshua Bengio introduced the architecture in “N-BEATS: Neural basis expansion analysis for interpretable time series forecasting” (2019, revised 2020). Start with the smallest unit.

A block takes a lookback window — the last L observations of a univariate series, typically some multiple of the forecast horizon H — and pushes it through a small stack of fully connected layers with ReLU activations. Out of that stack come two vectors of coefficients, not one:

input:  x ∈ ℝ^L        (the lookback window)

        h = FC_ReLU stack (x)

        θ_b = Linear_b(h)     coefficients for the BACKcast
        θ_f = Linear_f(h)     coefficients for the FOREcast

output: x̂ = B_b · θ_b   ∈ ℝ^L    (a reconstruction of the input window)
        ŷ = B_f · θ_f   ∈ ℝ^H    (the forecast)

The lookback length L is a hyperparameter and it is a more consequential one than it looks. Set it to two horizons and the model cannot see an annual cycle when forecasting a quarter ahead; set it to seven and each series yields far fewer training windows, because the first L observations of every series produce no target at all. The paper treats the lookback multiple as something to ensemble over rather than to tune, which is a tacit admission that no single value is right across a heterogeneous catalogue.

The block predicts the future and re-predicts the past it was given. That second output, the backcast, has no value to the user. It exists so the block can say what part of the input it believes it has accounted for, and that claim is what makes the stacking work.

The basis expansion

The B matrices above are the basis. Instead of the network emitting H numbers that are the forecast, it emits a small number of coefficients and the forecast is a linear combination of fixed or learned basis vectors. This is the same idea as writing a curve as a sum of polynomials or a sum of sinusoids: you constrain the space of shapes the output can take, and the network only chooses where in that space to sit.

In the generic configuration the basis is itself learned — B_f is just another linear layer with no structure imposed — and the paper is explicit that this configuration contains no time-series-specific components at all. It is fully connected layers and residual connections, and that is the headline claim: standard deep learning primitives, no seasonal differencing, no ARIMA order selection, no hand-built calendar features.

The reason the constraint helps rather than hurts is that the number of coefficients is much smaller than the horizon. A 48-step forecast built from 8 coefficients cannot wander; it can only produce shapes the basis spans. That is a strong regulariser applied at exactly the point where a forecaster tends to overfit, which is the far end of the horizon.

Doubly residual stacking

Blocks are chained, and the chaining is where the architecture’s name comes from. Each block passes on the part of the input it did not explain, and its forecast is added to a running total:

x_1 = x                        (first block sees the raw window)

for each block ℓ:
    x̂_ℓ, ŷ_ℓ = block_ℓ(x_ℓ)
    x_{ℓ+1}  = x_ℓ − x̂_ℓ        ← backward residual: what is left to explain
    Y        = Y + ŷ_ℓ           ← forward residual: forecasts accumulate

Two residual streams running in opposite directions, which is what doubly residual means. The backward stream is a sequential decomposition: block one removes what it can model, block two works on the remainder, and so on. The forward stream is an additive ensemble built during the forward pass.

The consequence worth holding on to is that the blocks are not interchangeable. Their order is meaningful, because each one sees a different signal — and if an early block over-explains the input, later blocks receive a residual with structure removed that was never really there, and their contribution is noise fitting. The paper stacks these many layers deep and relies on the residual path to keep the gradients usable, exactly as a ResNet does in vision.

The interpretable configuration

The second configuration replaces the learned basis with two fixed ones, arranged into two stacks:

• A trend stack whose basis is a small-degree polynomial in normalised time — 1, t, t², … up to a low maximum. Because the degree is small, the only shapes this stack can emit are slow, monotone-ish curves. It cannot produce a seasonal wiggle even if that would reduce the loss.
• A seasonality stack whose basis is a Fourier series — sines and cosines at harmonics of the seasonal period. Symmetrically, it can only emit periodic shapes.

Because the trend stack runs first and passes its residual on, the output decomposes: the trend stack’s accumulated forecast is the trend component, and the seasonality stack’s is the seasonal component. You get a decomposition of the same kind that STL produces, except that it falls out of the architecture rather than being computed beforehand.

The interpretability is real but limited in a specific way: it tells you how the forecast was composed, not why. There is nothing in it that attributes the trend to a driver, and the split between the two stacks depends on where the first stack’s polynomial degree was capped. Move that cap and the same series decomposes differently.

What it does not do

• It is univariate and takes no covariates. The input is a window of the series and nothing else — no price, no promotion, no holiday flag. This is a design decision, not an oversight, and it is the main reason to reach for the Temporal Fusion Transformer instead when you have known-future inputs. The follow-up work known as NBEATSx adds exogenous inputs; the original does not have them.
• It produces a point forecast. The architecture as published emits H numbers. Prediction intervals require either a quantile loss variant or an ensemble spread, and the ensemble spread is not a calibrated interval.
• The reported competition result is an ensemble. The M4 numbers come from ensembling across lookback lengths, loss functions and random initialisations rather than from one trained network. Check the paper’s experimental section for the exact ensemble size before budgeting for it — a single model is a materially different proposition to reproduce.
• It needs many series. Being a global model, it wants a catalogue, not a series. The observations-to-parameters arithmetic on when a classical model wins applies to it directly.