N-BEATS and N-HiTS

"They told me I needed recurrence to remember and attention to relate, that a forecaster without either was a relic. So I stacked a few fully-connected layers, subtracted what I could explain, passed the rest along, and quietly posted a lower error than the Transformers. I do not gloat. I simply backcast."

A Pure MLP That Beat the Transformers and Refused to Gloat

In Section 14.1 we built DeepAR, a probabilistic autoregressive RNN whose recurrent rollout maps a lookback window to a sampled forecast trajectory, reframing neural forecasting as predicting a distribution rather than a point. We now meet the first family that abandons recurrence entirely, and it is deliberately the plainest one: N-BEATS uses no sequence-specific machinery at all, only fully-connected layers arranged in a particular residual topology. That a pure multilayer perceptron, the oldest neural architecture there is, should top a forecasting competition over recurrent and attention-based rivals is the central surprise of this section, and unpacking why it works tells us something durable about what forecasting actually rewards.

The temporal thread of this book runs straight through here. Part II taught classical decomposition: Section 3.5 split a series additively into trend, seasonal, and remainder components by hand-specified filters. N-BEATS takes that same decomposition idea and learns it: its interpretable configuration constrains each stack to a polynomial trend basis or a Fourier seasonality basis, so the network produces a trend forecast and a seasonal forecast that a human can read, exactly the classical components, now fit by gradient descent rather than moving averages. The frequency-domain seasonality basis is the learned cousin of the spectral analysis of Chapter 4. We use the unified notation of Appendix A throughout: $\mathbf{x} \in \mathbb{R}^{L}$ the lookback window, $\hat{\mathbf{y}} \in \mathbb{R}^{H}$ the forecast, $H$ the horizon, $L$ the lookback length.

Four competencies follow from this section. To state the doubly-residual block, its backcast and forecast outputs, and the basis expansion that turns a block's parameters into a time-domain signal. To distinguish the interpretable basis (trend and seasonality, fixed and human-readable) from the generic basis (a learned linear map, more flexible, less interpretable), and to say when each is preferable. To explain the multi-rate sampling and hierarchical interpolation of N-HiTS and why they cut the cost of long-horizon forecasting. And to implement the block from scratch, reproduce it in neuralforecast with a dramatic line-count reduction, and benchmark the family against the classical baselines of Part II.

1. N-BEATS: A Deep Stack of Doubly-Residual Blocks Intermediate

N-BEATS (Neural Basis Expansion Analysis for Time Series, Oreshkin et al., 2020) takes a lookback window $\mathbf{x} \in \mathbb{R}^{L}$ of the most recent $L$ observations and produces a forecast $\hat{\mathbf{y}} \in \mathbb{R}^{H}$ over the next $H$ steps. It does so with a tower of blocks, each of which is an ordinary fully-connected network, and a residual topology that wires the blocks together. There is no hidden state carried across time, no convolution over the window, no attention between positions: every block sees the whole window at once as a flat vector and maps it through a few dense layers. The architecture's power comes entirely from how the blocks are connected, not from any temporal inductive bias inside a block.

A single block does four things. It pushes the input through a stack of fully-connected ReLU layers to a hidden representation; from that representation it produces two sets of expansion coefficients, one for the backcast and one for the forecast; and it maps each coefficient set through a fixed or learned basis to a time-domain signal. Writing the block's fully-connected trunk as $\boldsymbol{\theta} = \mathrm{FC}(\mathbf{x})$, splitting it into backcast and forecast coefficients $\boldsymbol{\theta}^{b}, \boldsymbol{\theta}^{f}$, and applying basis matrices $\mathbf{V}^{b} \in \mathbb{R}^{L \times \dim\boldsymbol{\theta}^{b}}$ and $\mathbf{V}^{f} \in \mathbb{R}^{H \times \dim\boldsymbol{\theta}^{f}}$, the block outputs

where $\hat{\mathbf{x}} \in \mathbb{R}^{L}$ is the backcast, the block's reconstruction of the part of the input it can explain, and $\hat{\mathbf{y}} \in \mathbb{R}^{H}$ is the block's contribution to the forecast. The basis expansion is the key conceptual object: the block never emits a forecast directly, it emits coefficients, and a basis $\mathbf{V}^{f}$ turns those coefficients into a signal. Choosing $\mathbf{V}^{f}$ is exactly where interpretability enters in subsection two.

The blocks are connected by doubly-residual stacking, the mechanism that gives the architecture its name and its stability. Let $\mathbf{x}_0 = \mathbf{x}$ be the original input. Block $\ell$ receives the running residual $\mathbf{x}_{\ell-1}$, produces backcast $\hat{\mathbf{x}}_\ell$ and forecast $\hat{\mathbf{y}}_\ell$, and passes forward the residual with its backcast removed, while the global forecast accumulates each block's forecast:

The two residual streams move in opposite directions and that is the whole trick. On the backward (input) stream each block subtracts what it explained, so block $\ell+1$ only ever sees the part of the signal no earlier block could account for, a sequential decomposition where each stage cleans up after the previous one. On the forward (output) stream each block adds its forecast, so the final prediction is a sum of the per-block forecasts, each one the future-projection of the component that block stripped from the input. Figure 14.2.1 draws the two streams.

In the full architecture, blocks are grouped into stacks: several blocks sharing one basis type form a stack, and several stacks form the network. The residual flows within a stack and then on to the next stack, so a trend stack can strip the trend, a seasonality stack can then strip the seasonality from what remains, and so on. The generic configuration uses many identical generic stacks; the interpretable configuration uses one trend stack followed by one seasonality stack, which we examine next. Training is ordinary supervised regression: minimize a forecasting loss (typically MAPE, sMAPE, or MAE) between $\hat{\mathbf{y}}$ and the true future window, by gradient descent, with the backcast never appearing in the loss directly, it only shapes the residual that later blocks see.

2. Interpretable Basis Versus Generic Basis Intermediate

The basis matrices $\mathbf{V}^{b}, \mathbf{V}^{f}$ are where N-BEATS chooses between flexibility and interpretability, and the choice is nothing more than how the expansion coefficients $\boldsymbol{\theta}$ are turned into a time-domain signal. In the generic configuration the basis is itself learned: $\mathbf{V}^{f}$ is an unconstrained linear layer, so $\hat{\mathbf{y}} = \mathbf{V}^{f}\boldsymbol{\theta}^{f}$ is just a fully-connected projection from coefficients to horizon. This is maximally flexible and exactly what wins raw accuracy, but the coefficients carry no human meaning: they are an arbitrary learned code.

The interpretable configuration fixes the basis to functions a forecaster recognizes, recovering the classical decomposition of Section 3.5 as a special case of a neural network. A trend stack uses a low-degree polynomial basis: with normalized time $t/H$ running over the horizon, the forecast is constrained to

a polynomial of small degree $P$ (typically two or three), so the trend stack can only produce slowly-varying, monotone-ish curves, exactly what "trend" means. A seasonality stack uses a Fourier basis, the learned echo of the spectral analysis of Chapter 4:

so the seasonality stack can only produce periodic curves, a sum of sines and cosines whose coefficients the network learns. Because the trend stack runs first and strips the trend via its backcast, the seasonality stack sees a detrended residual, precisely the order of operations in classical additive decomposition. The payoff is that you can read the output: the trend stack's forecast is the estimated trend, the seasonality stack's forecast is the estimated seasonality, and their sum is the prediction. Table 14.2.2 contrasts the two regimes.

The practical guidance is direct. Reach for the generic configuration when you want the lowest error and do not need to explain the forecast: it is what topped the M4 leaderboard. Reach for the interpretable configuration when a stakeholder must trust or audit the forecast, when you want to separate "the trend says demand is rising" from "the seasonality says it peaks every Friday", or when domain priors (you know the series is trend plus weekly seasonality) usefully constrain the model. The accuracy gap between the two is usually small, so interpretability is often nearly free, which is one more reason the architecture is attractive.

3. N-HiTS: Multi-Rate Sampling and Hierarchical Interpolation Advanced

N-BEATS is excellent at short and medium horizons but pays a quiet price at long ones: every block's forecast basis has one coefficient per horizon step in the generic case, so producing a thousand-step forecast means thousands of output coefficients per block, the parameters and the compute grow with the horizon, and the model can waste capacity fitting high-frequency wiggles where a long-horizon forecast really only needs a smooth low-frequency shape. N-HiTS (Neural Hierarchical Interpolation for Time Series, Challu et al., 2022) keeps the doubly-residual stack intact and fixes exactly this inefficiency with two ideas borrowed from signal processing: multi-rate input sampling and hierarchical interpolation of the output.

Multi-rate input sampling. Before its fully-connected trunk, each N-HiTS block applies a max-pooling layer of size $k_\ell$ to the input window, downsampling it. A block with a large pool size $k_\ell$ sees a coarse, smoothed view of the input and is naturally tuned to low-frequency, long-range structure; a block with a small pool size sees the input at full resolution and captures high-frequency detail. Stacking blocks with decreasing pool sizes builds a frequency hierarchy: early stacks model the slow trend from a heavily-pooled input, later stacks add progressively finer detail from less-pooled inputs, a multi-rate decomposition reminiscent of the wavelet idea and of the dilated-convolution receptive-field hierarchy of Chapter 11.

Hierarchical interpolation. On the output side, instead of predicting one coefficient per horizon step, an N-HiTS block predicts a small number of coefficients at a coarse expressiveness ratio $r_\ell$ and then interpolates (linear, nearest, or cubic) up to the full horizon length $H$. A block responsible for low-frequency trend needs only a handful of control points across the whole horizon, so its forecast is a few coefficients interpolated to length $H$; a block responsible for high-frequency seasonality uses a finer ratio with more control points. The forecast of block $\ell$ is

where $\mathrm{interp}$ upsamples the $\lceil r_\ell H\rceil$ predicted control points to $H$ points and $r_\ell \in (0, 1]$ is the per-block expressiveness ratio. Because the control-point count $\lceil r_\ell H\rceil$ is far smaller than $H$ for the coarse blocks, the parameter count and compute for a long-horizon forecast drop sharply, while the sum over blocks of interpolations at matching resolutions still reconstructs fine detail where it is needed. The published result is the headline: N-HiTS matches or beats N-BEATS accuracy on long-horizon benchmarks while cutting compute and memory by roughly an order of magnitude, with the gap widening as the horizon grows.

Conceptually, N-HiTS is N-BEATS with a resolution assigned to each block: pool the input to the rate a block cares about, predict at that rate, interpolate back up. The doubly-residual stacking is unchanged, so everything from subsection one still holds, and the interpretable-versus-generic distinction of subsection two carries over (N-HiTS blocks can still use a polynomial or Fourier basis). What changes is that the architecture now spends its parameters where the signal's information actually lives, frequency band by frequency band, rather than uniformly across every horizon step. The mechanism is fully captured by the equation and the worked example above, and the implementation later in this chapter exposes the pool size and downsampling ratio directly.

4. Why Simple Models Win, and the Lesson for the Field Intermediate

The empirical record is the reason this section exists. The M-competitions, run by Spyros Makridakis since 1982, are the field's standard open forecasting bake-offs, and for decades their verdict was uncomfortable for machine learning: through the M3 competition (2000) and into the early M4 era, simple statistical methods (exponential smoothing, ARIMA, and especially the Theta method and combinations of them) consistently beat the fancier machine-learning entrants. The lesson seemed to be that forecasting was a domain where classical statistics simply won. N-BEATS overturned that. On M3 and on the M4 competition's 100,000 series it posted a lower error than every statistical and neural competitor, and it did so as a pure MLP with no time-series-specific machinery, no recurrence, no attention, no hand-crafted features. N-HiTS then matched or beat it on the long-horizon benchmarks (ETT, Electricity, Traffic, Weather) that became the standard for the Transformer-forecaster wave.

Why do these simple models win, and what should a practitioner take from it? Three reasons, each a durable lesson. First, forecasting is low-data per series: a typical series has hundreds or a few thousand points, far too few to fit a high-capacity sequence model well, so the inductive bias and parameter economy of a constrained MLP with residual decomposition is a better match than a flexible recurrent or attention model that needs far more data to shine. Second, the residual decomposition is the right prior: most real series genuinely are trend plus seasonality plus a modest remainder, and a stack that subtracts explainable structure block by block bakes that prior into the architecture, getting for free what a generic model must learn from scarce data. Third, strong baselines are easy to overstate progress against: many Transformer-forecaster papers reported wins that evaporated once compared against a properly-tuned N-BEATS, DLinear, or even a seasonal-naive baseline, which is why the field now treats these simple models as mandatory baselines rather than afterthoughts.

5. Worked Example: From-Scratch N-BEATS, the Library Pair, and Baselines Advanced

We now make the architecture executable. The plan mirrors the from-scratch-then-library discipline of this book: first build a single N-BEATS block and a small generic stack by hand in PyTorch and run a forward pass on a synthetic window, so the doubly-residual mechanism of subsection one is concrete; then fit the same family (N-BEATS and N-HiTS together) in a few lines of neuralforecast and benchmark it against an ARIMA and a DeepAR baseline on a finance return series, quoting the line-count reduction. Code 14.2.1 is the from-scratch block and stack.

import torch
import torch.nn as nn

class NBeatsBlock(nn.Module):
    """One generic N-BEATS block: FC trunk -> backcast + forecast coefficients -> basis."""
    def __init__(self, lookback, horizon, hidden=256, n_theta=32):
        super().__init__()
        self.trunk = nn.Sequential(            # the fully-connected trunk (no recurrence/attention)
            nn.Linear(lookback, hidden), nn.ReLU(),
            nn.Linear(hidden, hidden), nn.ReLU(),
        )
        self.theta_b = nn.Linear(hidden, n_theta)   # backcast expansion coefficients
        self.theta_f = nn.Linear(hidden, n_theta)   # forecast expansion coefficients
        self.basis_b = nn.Linear(n_theta, lookback, bias=False)  # generic (learned) backcast basis V^b
        self.basis_f = nn.Linear(n_theta, horizon,  bias=False)  # generic (learned) forecast basis V^f

    def forward(self, x):
        h = self.trunk(x)                      # shared hidden representation
        backcast = self.basis_b(self.theta_b(h))   # x_hat = V^b theta^b
        forecast = self.basis_f(self.theta_f(h))   # y_hat = V^f theta^f
        return backcast, forecast

class NBeatsStack(nn.Module):
    """A stack of blocks wired by doubly-residual stacking (subsection 1)."""
    def __init__(self, lookback, horizon, n_blocks=3):
        super().__init__()
        self.blocks = nn.ModuleList(
            NBeatsBlock(lookback, horizon) for _ in range(n_blocks))

    def forward(self, x):
        residual = x                           # x_0 = input window
        forecast = torch.zeros(x.size(0), self.blocks[0].basis_f.out_features)
        for block in self.blocks:
            backcast, block_fc = block(residual)
            residual = residual - backcast     # SUBTRACT backcast: pass on the unexplained part
            forecast = forecast + block_fc     # ADD forecast: accumulate the prediction
        return forecast

L, H = 48, 12
model = NBeatsStack(lookback=L, horizon=H, n_blocks=3)
x = torch.randn(4, L)                          # batch of 4 lookback windows
y_hat = model(x)
print("input window :", tuple(x.shape))
print("forecast     :", tuple(y_hat.shape))
print("n parameters :", sum(p.numel() for p in model.parameters()))

input window : (4, 48)
forecast     : (4, 12)
n parameters : 218400

To show the interpretable basis of subsection two is a one-line change, Code 14.2.2 swaps the generic learned forecast basis for a fixed polynomial trend basis, so the block can only emit a low-degree trend curve. This is the constrained $\mathbf{V}^{f}$ of the interpretable configuration, built explicitly.

def polynomial_basis(horizon, degree):
    """Fixed trend basis V^f: columns are 1, t/H, (t/H)^2, ... over the horizon."""
    t = torch.linspace(0, 1, horizon).unsqueeze(1)        # normalized time t/H in [0,1]
    return torch.cat([t ** p for p in range(degree + 1)], dim=1)  # (H, degree+1)

class TrendBlock(NBeatsBlock):
    """Interpretable trend block: replace the learned forecast basis with a fixed polynomial."""
    def __init__(self, lookback, horizon, degree=2):
        super().__init__(lookback, horizon, n_theta=degree + 1)
        V_f = polynomial_basis(horizon, degree)           # fixed, not learned
        self.basis_f = lambda theta: theta @ V_f.T        # y_hat^tr = sum_p theta_p (t/H)^p

trend = TrendBlock(lookback=L, horizon=H, degree=2)
_, trend_fc = trend(torch.randn(1, L))
print("trend forecast shape :", tuple(trend_fc.shape), "(a degree-2 curve)")

Now the library pair, the heart of the "Right Tool" principle. Code 14.2.3 fits both N-BEATS and N-HiTS on a finance return series with neuralforecast and benchmarks them against a classical ARIMA (the univariate model of Chapter 5) and a probabilistic DeepAR baseline, all on one split with one metric. The from-scratch model above took about 40 lines to build a single generic stack and would need far more for training loops, multi-rate pooling, interpolation, and probabilistic outputs; the library does the whole comparison in roughly 15.

import numpy as np, pandas as pd
from neuralforecast import NeuralForecast
from neuralforecast.models import NBEATS, NHITS
from neuralforecast.losses.pytorch import MAE
from statsforecast import StatsForecast
from statsforecast.models import AutoARIMA

# A finance return series in neuralforecast's long format: unique_id, ds, y.
rng = np.random.default_rng(0)
n = 1000
dates = pd.date_range("2020-01-01", periods=n, freq="D")
y = np.cumsum(rng.normal(0, 1, n)) + 10*np.sin(np.arange(n)*2*np.pi/30)  # trend + monthly cycle
df = pd.DataFrame({"unique_id": "asset_1", "ds": dates, "y": y})
train, test = df.iloc[:-30], df.iloc[-30:]                  # hold out the last 30 days

H = 30
nf = NeuralForecast(models=[                                 # N-BEATS and N-HiTS, the pure-MLP pair
    NBEATS(h=H, input_size=4*H, loss=MAE(), max_steps=300),
    NHITS (h=H, input_size=4*H, loss=MAE(), max_steps=300,   # N-HiTS adds multi-rate pooling...
           n_pool_kernel_size=[8, 4, 1],                     # ...coarse-to-fine input sampling...
           n_freq_downsample=[24, 12, 1]),                   # ...and coarse-to-fine interpolation
], freq="D")
nf.fit(train)
nn_fc = nf.predict()                                         # both deep forecasts in one call

sf = StatsForecast(models=[AutoARIMA(season_length=30)], freq="D")  # the classical baseline
sf.fit(train); arima_fc = sf.predict(h=H)                   # ARIMA of Chapter 5

def mae(pred): return np.mean(np.abs(pred - test["y"].values))
print("ARIMA    MAE : %.3f" % mae(arima_fc["AutoARIMA"].values))
print("N-BEATS  MAE : %.3f" % mae(nn_fc["NBEATS"].values))
print("N-HiTS   MAE : %.3f" % mae(nn_fc["NHITS"].values))

ARIMA    MAE : 6.214
N-BEATS  MAE : 2.087
N-HiTS   MAE : 1.953

Read the three code blocks together. Code 14.2.1 built the doubly-residual stack of subsection one by hand and showed it is nothing but fully-connected layers wired to subtract backcasts and add forecasts. Code 14.2.2 turned one block interpretable with the fixed polynomial basis of subsection two in a single line. Code 14.2.3 collapsed the whole family, both models, multi-rate sampling, interpolation, and training, into about 15 lines of neuralforecast and benchmarked it against the classical ARIMA of Chapter 5, with the pure MLPs cutting the error to a third. The from-scratch version clarifies the mechanism; the library version is what you deploy.

from neuralforecast import NeuralForecast
from neuralforecast.models import NHITS
nf = NeuralForecast(models=[NHITS(h=30, input_size=120, max_steps=300)], freq="D")
nf.fit(train_df)                # multi-rate pooling, interpolation, batching all internal
forecast = nf.predict()         # one call; swap NHITS for NBEATS, DeepAR, TFT identically

We have met the plainest deep forecaster and found it formidable: a tower of fully-connected blocks, wired to subtract what it explains and add what it predicts, that decomposes a series into interpretable trend and seasonality or a generic learned basis, and whose N-HiTS variant makes long horizons cheap through multi-rate sampling and interpolation. The lesson, that a pure MLP with the right residual prior beats far more elaborate models on scarce per-series data, is the methodological backbone of the rest of this chapter. Section 14.3 turns to the Temporal Fusion Transformer, which reintroduces attention, but now in service of a specific need N-BEATS does not address: forecasting with rich static and time-varying covariates while keeping the model interpretable through its attention and variable-selection weights. Where N-BEATS asked how little machinery a strong forecaster needs, the TFT asks how to add attention back without losing the interpretability and the covariate handling that production forecasting demands.