Machine Learning for Time Series

Contents

• Introduction and overview
• 0. Review: probability, least squares and Fourier analysis
• 1. Representing a time series: stationarity, spectra and scales
• 2. Linear models, state space and forecasting
• 3. Extracting and recognising patterns
• 4. Sparse representations and dictionary learning
• 5. Detecting changes and anomalies
• 6. Repairing data: denoising and missing values
• 7. Extracting and selecting relevant features
• 8. Multivariate series and graphs
• 9. Neural networks for time series
• 10. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student with an M1 in mathematics, computer science or physics, who has met probability and a first course in machine learning, and who now faces data indexed by time — sensors, physiology, energy, economics.

Empirical first: on the synthetic Cylinder–Bell–Funnel benchmark, a 1-nearest-neighbour classifier with the Euclidean distance misclassifies about 13 % of the test series, and the same classifier with dynamic time warping misclassifies none of them (code of ch. 3). The Euclidean distance compares sample t with sample t, while the patterns that define the classes are shifted and stretched in time.

That small experiment is the whole course in miniature. A time series is not a vector of independent coordinates: neighbouring samples are dependent, patterns move, regimes change, values go missing, and several channels interact. Each chapter chooses a representation that takes the temporal structure into account — autocorrelations and spectra, linear dynamics, warping paths, sparse codes, segmentations, features, graphs — and turns it into an algorithm.

The page leans on the optimisation of Convex Optimization and Applications in Machine Learning and the statistics of Introduction to Statistical Learning. Its graph chapter opens the door to Graphs in Machine Learning, and its detection chapter to Detection Theory and Industrial Applications.

What you will be able to do

• Compute and interpret the autocorrelation, the periodogram, the spectrogram and the wavelet coefficients of a series.
• Fit autoregressive and state-space models, forecast with them, and evaluate any forecaster without leaking the future.
• Compare series with dynamic time warping, bound it from below, and build nearest-neighbour and shapelet classifiers.
• Discover repeated motifs and isolated discords with the matrix profile.
• Learn sparse and convolutional dictionaries that extract recurring patterns from long signals.
• Segment a signal by penalised change-point detection (optimal partitioning, PELT, kernel costs) and choose the penalty.
• Design and evaluate detectors for point, contextual and collective anomalies.
• Denoise by Wiener filtering and wavelet thresholding, and impute missing values by splines, Kalman smoothing and low-rank completion.
• Extract interpretable features from series and select them without selection bias.
• Represent a multivariate series on a graph, filter and interpolate it, and learn the graph from the data.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com starts from applications (health, economics, advertising…) in which the collected data take the form of time series. The fundamental issue, it says, is to choose a representation that accounts as well as possible for the temporal information.

The course aims to give an overview of the machine-learning techniques used to study time series, for tasks such as pattern extraction and recognition, anomaly detection, prediction and interpolation. It is largely illustrated on real data and problems from current challenges, and insists on understanding and interpreting the data.

master-mva.com announces 6 lecture sessions of 3 h, 3 lab sessions of 3 h and the defences, and states that the course does not accept external auditors. The Université Paris-Saclay course detail adds only the nature of the assessment. The table below maps every official item to a section of this page.

Beyond the syllabus, this page adds the Kalman smoother as a unifying tool for forecasting and imputation (ch. 2 and 6), neural networks for time series (ch. 9, self-study) and a chapter of research directions (ch. 10, self-study).

Where this course fits in the programme

Official prerequisites (master-mva.com): probability and statistics, and basic notions of machine learning and signal processing.

• Refreshers — Probability and Statistics (Gaussian vectors, conditional expectation, estimation), Linear Algebra, Analysis and Optimization (least squares, eigenvalues, gradient methods) and Scientific Python, NumPy and PyTorch (every code block of this page).
• Related courses of the same semester — Convex Optimization (the lasso and proximal methods of ch. 4 and 7), Introduction to Statistical Learning, Introduction to Probabilistic Graphical Models and Deep Generative Models (state-space models, Gaussian graphical models), Computational Statistics, Statistical Learning with Extreme Values (rare events and anomalies), Deep Learning and Signal Processing and Mathematical Methods for Neuroscience.
• Courses that build on it in semester 2 — Graphs in Machine Learning, Detection Theory and Industrial Applications, Sequential Learning, Kernel Methods for Machine Learning (kernel change-point costs), Bayesian Machine Learning, Audio Signal Processing: Time–Frequency Analysis, Functional Brain Imaging and Brain–Computer Interfaces, Remote Sensing Data and Biostatistics.

Course outline

• Review: probability, least squares and Fourier analysis — Gaussian conditioning, ridge regression, the discrete Fourier transform and the convolution theorem.
• Representing a time series — stationarity, autocorrelation, spectral density, periodogram, spectrogram, Haar wavelets, trend and seasonality.
• Linear models, state space and forecasting — AR and ARMA models, Yule–Walker equations, best linear prediction, the Kalman filter, forecasting as supervised learning and its honest evaluation.
• Extracting and recognising patterns — dynamic time warping and its lower bounds, nearest-neighbour classification, the matrix profile, shapelets.
• Sparse representations and dictionary learning — OMP, ISTA and FISTA, recovery guarantees, MOD and K-SVD, convolutional sparse coding.
• Detecting changes and anomalies — cost functions, optimal partitioning, PELT, kernel costs, penalties; anomaly types, residual, isolation and sequential detectors.
• Repairing data — Wiener filtering, wavelet thresholding, missingness mechanisms, splines, Kalman smoothing, low-rank completion.
• Extracting and selecting relevant features — interpretable features, random convolutional kernels, filter, wrapper and embedded selection.
• Multivariate series and graphs — graph Fourier transform, graph filters, interpolation on a graph, learning the graph.
• Neural networks for time series — dilated convolutions, recurrent networks and Transformers, and the baselines they must beat (self-study).
• Research directions and open problems — foundation models, evaluation, irregular sampling, spatio-temporal graphs (self-study).

Session plan

36 h in class = 12 sessions of 3 h. The volume is that of the Université Paris-Saclay programme table (lectures 27 h, labs 9 h). master-mva.com organises the course as 6 lectures of 3 h, 3 labs of 3 h and the defences; the plan below spreads the 27 h of lectures over 9 sessions and keeps the 3 labs. The split into sessions is a proposal.

Level 1 = foundations and first examples, 2 = core models and algorithms, 3 = advanced material close to research; the level never decreases.

Self-study, outside the 36 h (about 6 h in total, plus the mini-project): ch. 9 Neural networks for time series (≈ 3 h, with Ex. 15) · ch. 10 Research directions and open problems (≈ 3 h, with Ex. 16).

How to study this module

• Read ch. 0 before session 1 and check that Proposition 0.2 (Gaussian conditioning) and Theorem 0.5 (convolution theorem) feel familiar: ch. 2, 4, 6 and 8 use them without warning.
• After each session, do the exercises listed in the Practice column, then rerun the code blocks of the chapter and change one parameter at a time.
• Keep one notebook per lab; the lab reports count for 25 % of the official grade, so write the observations as you go.
• Choose the mini-project early (session 8 at the latest): most subjects reuse a lab.
• Python environment: NumPy, SciPy, scikit-learn, NetworkX and PyTorch (CPU) suffice for every block; all data are synthetic, with fixed seeds.

0. Review: probability, least squares and Fourier analysis

Three M1 tools carry most of the course: conditioning of Gaussian vectors (prediction, Kalman filtering, imputation), least squares (autoregressions, dictionaries, features) and the discrete Fourier transform (spectra, convolutions, filters). This chapter states them in the form used later, with short proofs.

The same orthogonality principle — the error is orthogonal to everything the predictor may use — gives the Yule–Walker equations (ch. 2), the Kalman filter (ch. 2) and the Wiener filter (ch. 6).

A linear filtering, a cross-correlation or a sliding inner product therefore costs \mathcal{O}(n \log n) instead of \mathcal{O}(n^2). This speeds up the matrix profile (ch. 3) and convolutional sparse coding (ch. 4).

1. Representing a time series: stationarity, spectra and scales

Before learning anything from a series, one must decide what to look at: the raw samples, their correlations across lags, their frequencies, or their variations at several scales. This chapter defines these representations, says when they are meaningful (stationarity) and how to estimate them.

Time series as data

Stationarity and autocorrelation

Spectral representation

The autocovariance describes dependence across lags; its Fourier transform describes how the variance is spread over frequencies.

White noise has a flat spectrum f = \sigma^2 / (2\pi). The AR(1) of Example 1.6 has f(\omega) = \frac{\sigma^2}{2\pi\,(1 - 2\phi\cos\omega + \phi^2)}, which concentrates the variance at low frequencies when \phi > 0. It follows by summing the two geometric series of \sum_h \phi^{|h|} e^{-ih\omega}.

The periodogram must therefore be smoothed. Welch’s method cuts the series into K overlapping windowed segments and averages their periodograms. The variance is divided by about K, and the frequency resolution, now set by the segment length, is coarser by the same factor.

import numpy as np
from scipy.signal import welch

rng = np.random.default_rng(0)
n, phi = 4096, 0.7
eps = rng.standard_normal(n)
ar = np.zeros(n)
for t in range(1, n):                       # AR(1): x_t = 0.7 x_{t-1} + eps_t
    ar[t] = phi * ar[t - 1] + eps[t]
x = ar + 1.5 * np.sin(2 * np.pi * 0.1 * np.arange(n))   # + a period-10 cycle

def acf(x, max_lag):                        # sample ACF (Definition 1.5)
    xc = x - x.mean()
    g = np.array([xc[: len(x) - h] @ xc[h:] for h in range(max_lag + 1)]) / len(x)
    return g / g[0]

print("sample ACF, lags 0-3:", acf(ar, 3).round(2))
print("theory 0.7^h        :", (phi ** np.arange(4)).round(2))
f, P = welch(x, fs=1.0, nperseg=250)        # averaged periodogram
print("frequency of the spectral peak:", f[np.argmax(P)], "cycle/sample")

The output gives the ACF (1, 0.70, 0.49, 0.35) against the theoretical (1, 0.70, 0.49, 0.34), and a spectral peak at f = 0.1 \text{ cycle per sample}, the hidden period of 10 samples, invisible on the raw ACF of the noisy mixture.

Time–frequency and multiscale representations

A spectrum averages over the whole series. Speech, gait or an electrocardiogram change over time, so one needs representations that are localised both in time and in frequency.

Wavelets resolve the trade-off differently: short windows for high frequencies, long ones for low frequencies. The simplest wavelet basis, Haar’s, is enough to state the main facts.

A piecewise-constant signal has nonzero Haar details only near its jumps: few large coefficients carry it. This sparsity is what wavelet denoising exploits in ch. 6.

Trend and seasonality

Many series are a slow trend plus a periodic seasonal component plus a stationary remainder, x_t = m_t + s_t + r_t. Differencing removes the first two.

Differencing thus removes trends but amplifies high frequencies: it should be used once or twice, not by reflex.

2. Linear models, state space and forecasting

To forecast a series, or to fill its gaps, one needs a model of how the present depends on the past. This chapter builds the linear one: autoregressions, their optimal predictors, the state-space form with the Kalman filter, and then forecasting as ordinary supervised learning with an honest evaluation.

Autoregressive and moving-average models

Replacing \gamma by the sample autocovariance (Definition 1.5) gives the Yule–Walker estimator \hat{\boldsymbol\phi} = \hat\Gamma_p^{-1} \hat{\boldsymbol\gamma}_p .

The 1/n normalisation keeps \hat\Gamma_p positive definite, and the fitted model is always causal. The estimator is consistent and asymptotically Gaussian; least squares on the lagged values is asymptotically equivalent.

State-space models and the Kalman filter

The state-space form writes every linear model — AR, ARMA, trend plus season, several sensors — as a hidden Markov chain observed with noise. One algorithm then forecasts, filters, smooths and handles missing values (ch. 6).

The innovations y_t - H m_{t \mid t-1} are independent \mathcal{N}(0, S_t) variables, so the likelihood of the whole series is a product of n Gaussian densities; maximising it fits F, Q and R. The Rauch–Tung–Striebel smoother adds a backward pass that conditions on the whole series: m_{t \mid n} = m_{t \mid t} + J_t\,(m_{t+1 \mid n} - m_{t+1 \mid t}), \qquad J_t = P_{t \mid t} F^\top P_{t+1 \mid t}^{-1}.

It follows from Proposition 0.2 applied to (s_t, s_{t+1}) given y_{1:t}, and from the Markov property: given s_{t+1}, the state s_t is independent of y_{t+1:n}. Chapter 6 uses it to impute gaps.

Forecasting as supervised learning

import numpy as np
from scipy.linalg import toeplitz, solve
from sklearn.linear_model import Ridge

rng = np.random.default_rng(1)
n = 2000
x = np.zeros(n)
for t in range(2, n):                       # AR(2) with phi = (0.6, -0.3)
    x[t] = 0.6 * x[t - 1] - 0.3 * x[t - 2] + rng.standard_normal()

def yule_walker(x, p):                      # Theorem 2.3 with sample autocovariances
    xc = x - x.mean()
    g = np.array([xc[: len(x) - h] @ xc[h:] for h in range(p + 1)]) / len(x)
    phi = solve(toeplitz(g[:p]), g[1:])
    return phi, g[0] - phi @ g[1:]

phi_hat, s2_hat = yule_walker(x[:1000], 2)
print("Yule-Walker on the first 1000 points:", phi_hat.round(3), round(s2_hat, 3))

L = 5                                       # lag embedding: row t holds x_{t-1}, ..., x_{t-5}
X = np.column_stack([x[L - k - 1 : n - k - 1] for k in range(L)])
y = x[L:]
mse = {"naive (last value)": [], "AR(2) Yule-Walker": [], "ridge on 5 lags": []}
for origin in range(1000, len(y), 100):     # rolling origin (Definition 2.10)
    te = slice(origin, origin + 100)
    phi, _ = yule_walker(x[: origin + L], 2)            # refit on the past only
    ridge = Ridge(alpha=1.0).fit(X[:origin], y[:origin])
    mse["naive (last value)"].append(np.mean((y[te] - X[te, 0]) ** 2))
    mse["AR(2) Yule-Walker"].append(np.mean((y[te] - X[te, :2] @ phi) ** 2))
    mse["ridge on 5 lags"].append(np.mean((y[te] - ridge.predict(X[te])) ** 2))
for name, v in mse.items():
    print(f"{name:20s} one-step MSE = {np.mean(v):.3f}")

The estimate is \hat{\boldsymbol\phi} = (0.578, -0.297) for the true (0.6, -0.3). The one-step errors are 1.567 for the naive forecast, 1.055 for the AR(2) model and 1.060 for the ridge regression on 5 lags.

Both models come close to the floor \sigma^2 = 1 of Theorem 2.4. The remaining gap of about 0.05 is within the sampling error of an average over about 1000 test points, whose standard deviation is close to \sqrt{2/1000} \approx 0.045. On a linear process, a flexible regressor cannot beat the right linear model; it can only match it.

3. Extracting and recognising patterns

Two heartbeats, two steps of a walker or two spoken words are “the same pattern” even when one is shifted, slower or locally stretched. This chapter builds distances that tolerate such deformations, uses them to classify series, and then searches long series for repeated patterns (motifs) and for short discriminative ones (shapelets).

Comparing two series

Metric indexing structures therefore do not apply directly; fast search relies on lower bounds instead.

In a nearest-neighbour search, a candidate whose lower bound exceeds the best distance found so far is discarded without computing DTW. The bound is presented in Keogh and Ratanamahatana (2005); cascades of such bounds let DTW search scale to trillions of subsequences (Rakthanmanon and coauthors, 2012).

Nearest-neighbour classification

DTW is not a metric, so the theorem is only a guide. Empirically, 1-NN with a DTW whose band is tuned by cross-validation was long the baseline to beat on the UCR archive of time-series classification data sets (Dau and coauthors, 2019). The large comparison of Bagnall and coauthors (2017) found that many published methods were not significantly more accurate than this baseline.

The Cylinder–Bell–Funnel data set (Saito, 1994) is synthetic: a plateau, a ramp up or a ramp down, of random start, length and height, in Gaussian noise. The code generates it and compares the two distances.

import numpy as np

def cbf(n_per_class, rng, length=128):      # Cylinder-Bell-Funnel (Saito, 1994)
    t, X, y = np.arange(length), [], []
    for label in range(3):
        for _ in range(n_per_class):
            a = rng.integers(16, 33)
            b = a + rng.integers(32, 97)
            chi = ((t >= a) & (t <= b)).astype(float)
            shape = [chi, chi * (t - a) / (b - a), chi * (b - t) / (b - a)][label]
            X.append((6 + rng.standard_normal()) * shape + rng.standard_normal(length))
            y.append(label)
    return np.array(X), np.array(y)

def dtw_to_all(q, C, w):
    """DTW between q and every row of C (squared costs, band half-width w)."""
    n, N = len(q), len(C)
    D = np.full((n + 1, n + 1, N), np.inf)
    D[0, 0] = 0.0
    for i in range(1, n + 1):                # recursion of Theorem 3.3
        for j in range(max(1, i - w), min(n, i + w) + 1):
            best = np.minimum(np.minimum(D[i - 1, j], D[i, j - 1]), D[i - 1, j - 1])
            D[i, j] = (q[i - 1] - C[:, j - 1]) ** 2 + best
    return np.sqrt(D[n, n])

rng = np.random.default_rng(2)
Xtr, ytr = cbf(10, rng)                      # 30 training series
Xte, yte = cbf(30, rng)                      # 90 test series
z = lambda X: (X - X.mean(1, keepdims=True)) / X.std(1, keepdims=True)
Xtr, Xte = z(Xtr), z(Xte)

d_euc = ((Xte[:, None, :] - Xtr[None, :, :]) ** 2).sum(-1)
d_dtw = np.array([dtw_to_all(q, Xtr, w=12) for q in Xte])
for name, d in [("Euclidean", d_euc), ("DTW, 10% band", d_dtw)]:
    print(f"1-NN {name:14s} error = {np.mean(ytr[d.argmin(1)] != yte):.3f}")

With 30 training series, the Euclidean 1-NN errs on 13.3 % of the test series and the DTW 1-NN on none: \text{error}_{\text{Euclid}} = 0.133, \qquad \text{error}_{\text{DTW}} = 0.000 .

Motifs and discords: the matrix profile

In a single long series, one asks which pattern repeats (a motif) and which stretch resembles nothing else (a discord). The matrix profile answers both questions at once.

All distances therefore follow from the sliding dot products \langle T_{i,m}, T_{j,m} \rangle and from sliding means and standard deviations (cumulative sums). For a fixed i, the dot products with all j form one cross-correlation, computed by FFT (Theorem 0.5). Along a diagonal they update in constant time: \langle T_{i+1,m}, T_{j+1,m} \rangle = \langle T_{i,m}, T_{j,m} \rangle - T_i T_j + T_{i+m} T_{j+m}.

This gives the exact matrix profile in \mathcal{O}(n^2) time and \mathcal{O}(n) memory (Yeh and coauthors, 2016; Zhu and coauthors, 2016). The block below uses the \mathcal{O}(n^2)-memory version, which is enough for 2000 samples.

import numpy as np

rng = np.random.default_rng(3)
n, m = 2000, 50
T = np.cumsum(rng.standard_normal(n)) * 0.1 + 0.3 * rng.standard_normal(n)
motif = 2 * np.sin(np.linspace(0, 3 * np.pi, m))
for pos in (300, 1400):                      # plant the same pattern twice
    T[pos : pos + m] += motif

S = np.lib.stride_tricks.sliding_window_view(T, m)        # all subsequences
Z = (S - S.mean(1, keepdims=True)) / S.std(1, keepdims=True)
corr = Z @ Z.T / m                                         # Pearson correlations
dist = np.sqrt(np.maximum(2 * m * (1 - corr), 0))          # Proposition 3.10
idx = np.arange(len(Z))
dist[np.abs(idx[:, None] - idx[None, :]) < m // 2] = np.inf   # exclusion zone
P = dist.min(1)                                            # the matrix profile
i = int(np.argmin(P))
j = int(np.argmin(dist[i]))
print("motif pair (start indices):", sorted((i, j)), " distance:", round(P[i], 2))
print("discord (most isolated subsequence) starts at:", int(np.argmax(P)))

The top motif is the pair of subsequences starting at 302 and 1402, two samples after the planted copies (the z-normalised pattern matches best with a small offset). The discord, at index 63, is simply the most unusual stretch of the random walk: on a series without a planted anomaly, the top discord is not an anomaly.

Shapelets

A class can be defined by a short pattern that appears somewhere in the series. Shapelets are such patterns, cut out of the training data.

With N series of length n, there are about N n^2 / 2 candidates, each compared with N series at a cost of order n^2. Exhaustive search is thus of order N^2 n^4 operations, which is why practical versions subsample the candidates or learn the shapelets by gradient descent (Grabocka and coauthors, 2014). The original shapelets of Ye and Keogh (2009) were used inside a decision tree.

The reward is interpretability: a shapelet is a piece of real data that one can plot and show to a domain expert.

import numpy as np
from sklearn.linear_model import LogisticRegression

rng = np.random.default_rng(4)
tri = np.r_[np.linspace(0, 2, 10), np.linspace(2, 0, 10)]   # class-1 pattern
sq = np.r_[np.zeros(3), 2 * np.ones(14), np.zeros(3)]      # class-0 pattern

def make(n, length=100):                    # pattern at a random position + noise
    X, y = 0.5 * rng.standard_normal((n, length)), rng.integers(0, 2, n)
    for k in range(n):
        p = rng.integers(0, length - 20)
        X[k, p : p + 20] += tri if y[k] else sq
    return X, y

def sdist(s, X):                            # shapelet distance (Definition 3.11)
    W = np.lib.stride_tricks.sliding_window_view(X, len(s), axis=1)
    return ((W - s) ** 2).mean(-1).min(-1)

def quality(d, y):                          # separation of the two classes
    a, b = d[y == 0], d[y == 1]
    return (a.mean() - b.mean()) ** 2 / (a.var() + b.var())

Xtr, ytr = make(40)
Xte, yte = make(300)
cands = [x[p : p + 20] for x in Xtr[:10] for p in range(0, 81, 2)]
scores = [quality(sdist(s, Xtr), ytr) for s in cands]
best = [cands[i] for i in np.argsort(scores)[::-1][:5]]   # keep the 5 best
F = lambda X: np.column_stack([sdist(s, X) for s in best]) # shapelet transform
clf = LogisticRegression().fit(F(Xtr), ytr)
print("shapelet transform + logistic regression, accuracy:", round(clf.score(F(Xte), yte), 3))
d = ((Xte[:, None, :] - Xtr[None, :, :]) ** 2).sum(-1)
print("Euclidean 1-NN on the raw series, accuracy:      ", round(np.mean(ytr[d.argmin(1)] == yte), 3))

The patterns sit at random positions, so the lock-step 1-NN reaches only 0.67, while five shapelets learned from 40 series give \text{accuracy} = 0.917 .

4. Sparse representations and dictionary learning

A recurring pattern — a heartbeat, a step, a spike — can be seen as an atom that the signal uses a few times. This chapter represents signals as sparse combinations of atoms, learns the atoms from the data, and then makes the representation shift-invariant with convolutions, which turns it into a pattern-extraction tool.

Sparse coding

For the union of the identity and the orthonormal Fourier basis in dimension n, \mu = 1/\sqrt n. With n = 64 the theorem guarantees recovery of any combination of at most 4 spikes and sinusoids (Exercise 10).

Dictionary learning

Fixed bases (Fourier, wavelets) are universal but not adapted to a given signal. Dictionary learning fits the atoms to the data (Tošić and Frossard, 2011, in the official references).

For millions of windows, the online algorithm of Mairal and coauthors (2010) updates the dictionary from mini-batches.

Convolutional sparse coding

Cutting a long series into windows creates an artefact: a pattern that straddles two windows is coded twice, at arbitrary shifts. The convolutional model codes the whole series at once, with atoms that may appear at any time.

One FISTA or ADMM iteration therefore costs a few FFTs, \mathcal{O}(K n \log n), instead of \mathcal{O}(K n L) for direct convolutions. The block below codes a synthetic signal built from two atoms with FISTA.

import numpy as np
from scipy.signal import find_peaks

rng = np.random.default_rng(5)
n, L, K = 1024, 32, 2
t = np.linspace(-1, 1, L)
atoms = np.stack([np.exp(-8 * t**2), np.sin(3 * np.pi * t) * np.exp(-3 * t**2)])
atoms /= np.linalg.norm(atoms, axis=1, keepdims=True)      # unit-norm atoms
Z = np.zeros((K, n))
for k in range(K):                                         # 5 activations per atom
    Z[k, rng.choice(n - L, 5, replace=False)] = rng.uniform(1, 3, 5)

Dh = np.fft.rfft(atoms, n)                                 # atoms in the Fourier domain
synth = lambda Z: np.fft.irfft((Dh * np.fft.rfft(Z, n)).sum(0), n)   # sum_k d_k * z_k
adj = lambda r: np.fft.irfft(np.conj(Dh) * np.fft.rfft(r, n), n)     # its adjoint
x = synth(Z) + 0.05 * rng.standard_normal(n)

lam, Lip = 0.2, (np.abs(Dh) ** 2).sum(0).max()            # Lipschitz constant (Prop. 4.9)
z = v = np.zeros((K, n))
tk = 1.0
for it in range(300):                                      # FISTA iterations
    u = v - adj(synth(v) - x) / Lip                        # gradient step
    z_new = np.sign(u) * np.maximum(np.abs(u) - lam / Lip, 0)   # soft threshold
    t_new = (1 + np.sqrt(1 + 4 * tk**2)) / 2
    v = z_new + (tk - 1) / t_new * (z_new - z)
    z, tk = z_new, t_new
for k in range(K):
    found, _ = find_peaks(z[k], height=0.5)
    print(f"atom {k}: true {np.flatnonzero(Z[k])}, found {found}")
print("relative residual:", round(np.linalg.norm(x - synth(z)) / np.linalg.norm(x), 3))

The ten activation times are recovered exactly, for both atoms. The relative residual of 0.23 reflects the shrinkage of the amplitudes by the \ell_1 penalty; refitting the amplitudes by least squares on the detected support removes this bias.

5. Detecting changes and anomalies

A gait recording switches from walking to running, a machine drifts before failing, a sensor freezes for an hour. The first two are changes: the law of the series is different after some times. The last one is an anomaly: a value or a stretch that does not conform to normal behaviour. This chapter treats change-point detection as an optimisation problem, following the official review of Truong, Oudre and Vayatis (2020), and then the detection of anomalies.

Change-point detection as penalised optimisation

Search methods

For c_{L_2} one can take \kappa = 0, since splitting never increases the cost (Proposition 5.3). The result is exact: PELT returns the same segmentation as optimal partitioning. When the number of changes grows linearly with n, its expected cost is \mathcal{O}(n) under mild conditions (Killick and coauthors, 2012).

The kernel cost is therefore the quadratic cost of the data mapped into \mathcal{H}. With a characteristic kernel, such as the Gaussian kernel, mean embeddings differ whenever distributions differ: the cost detects changes in variance or shape, not only in the mean (Harchaoui and Cappé, 2007; Arlot, Celisse and Harchaoui, 2019).

import numpy as np

rng = np.random.default_rng(6)
means, lengths = [0, 2, 0.5, 3, 1], [200, 150, 250, 100, 300]
y = np.concatenate([m + rng.standard_normal(l) for m, l in zip(means, lengths)])
n = len(y)
S1, S2 = np.r_[0, np.cumsum(y)], np.r_[0, np.cumsum(y**2)]
cost = lambda a, b: S2[b] - S2[a] - (S1[b] - S1[a]) ** 2 / (b - a)  # c_L2(y[a:b])

def pelt(beta, min_size=5, prune=True):
    F = np.full(n + 1, np.inf)
    F[0] = -beta
    last, R, evals = np.zeros(n + 1, int), [0], 0
    for t in range(min_size, n + 1):
        ok = [s for s in R if t - s >= min_size]
        vals = np.array([F[s] + cost(s, t) + beta for s in ok])
        evals += len(ok)
        F[t], last[t] = vals.min(), ok[int(vals.argmin())]
        if prune:                            # Theorem 5.6 with kappa = 0
            R = [s for s in R if t - s < min_size] + \
                [s for s, v in zip(ok, vals) if v - beta <= F[t]]
        R.append(t)
    cps, t = [], n
    while t > 0:                             # backtrack the optimal segmentation
        t = last[t]
        cps.append(t)
    return [int(c) for c in sorted(cps)[1:]], evals

beta = 2 * np.log(n)                         # BIC-type penalty, sigma = 1
print("true change points:", np.cumsum(lengths)[:-1].tolist())
for prune in (False, True):
    cps, evals = pelt(beta, prune=prune)
    print(f"{'PELT' if prune else 'optimal partitioning':20s} {cps}  cost evaluations: {evals}")

Both searches return the change points (201, 350, 600, 700) for the true (200, 350, 600, 700). Optimal partitioning needs 492 532 cost evaluations, PELT 100 618: the same exact answer, about five times cheaper here, and the gap grows with n.

Anomaly detection

The main families reuse the previous chapters. Prediction-based detectors score the residual of a forecaster (ch. 2). Distance-based ones score the distance to the nearest neighbour, such as the matrix-profile discords of ch. 3. Reconstruction-based ones score the error of a dictionary or a low-rank model (ch. 4).

With k = 3 the per-sample rate is 0.0027, that is 27 false alarms expected in 10^4 normal samples. To keep a 5 % chance of any false alarm over 10^4 samples, the threshold must rise to k \approx 4.56 .

For a Gaussian mean shift from \mu_0 to \mu_1, \ell_t is proportional to y_t - (\mu_0 + \mu_1)/2. The threshold h trades the mean time between false alarms against the detection delay; Lorden (1971) proved that CUSUM is asymptotically optimal for this trade-off. The theory is developed in Detection Theory and Industrial Applications.

import numpy as np
from sklearn.ensemble import IsolationForest

rng = np.random.default_rng(7)
n, p = 3000, 50                              # period-50 cycle + AR(1) noise
noise = np.zeros(n)
for t in range(1, n):
    noise[t] = 0.5 * noise[t - 1] + 0.3 * rng.standard_normal()
x = np.sin(2 * np.pi * np.arange(n) / p) + noise
x[700] += 3.0                                # point anomaly
x[1500:1525] = x[1475:1500]                  # contextual: half a cycle repeated
x[2300:2400] = np.sin(4 * np.pi * np.arange(100) / p) + noise[2300:2400]  # collective

# (a) residuals of a seasonal-naive forecast, standardised by a robust scale
r = x[p:] - x[:-p]
z = np.abs(r - np.median(r)) / (1.4826 * np.median(np.abs(r - np.median(r))))
def regions(idx, gap=10):                    # group flagged times into intervals
    cuts = np.flatnonzero(np.diff(idx) > gap)
    return [(int(a), int(b)) for a, b in zip(idx[np.r_[0, cuts + 1]], idx[np.r_[cuts, -1]])]
print("(a) residual |z| > 4 in:", regions(np.flatnonzero(z > 4) + p))

# (b) isolation forest on sliding windows of one period
W = np.lib.stride_tricks.sliding_window_view(x, p)[::5]
s = -IsolationForest(random_state=0).fit(W).score_samples(W)
print("(b) isolation forest, 5 most anomalous windows start at:", np.sort(np.argsort(s)[-5:] * 5))

The residual detector (a) flags the spike at 700, the contextual anomaly around 1508–1519 and parts of the collective one (2315–2335). It also flags 750, 1560–1567 and 2418–2433: these are echoes, because the seasonal forecast \hat x_t = x_{t-50} reuses the anomalous values one period later. The isolation forest (b) flags only windows that overlap the collective anomaly (starts 2285 to 2345).

Each detector sees the anomalies its representation makes visible; combining a residual score with a window score covers the three types.

6. Repairing data: denoising and missing values

Real recordings are noisy and incomplete: sensors add measurement noise, connections drop, batteries die. This chapter removes noise by filtering and by thresholding in a sparse basis, then fills missing values by interpolation, by the Kalman smoother of ch. 2, and by low-rank completion when several channels are recorded together.

Denoising

Throughout, y_t = x_t + \varepsilon_t: a signal of interest plus noise, and the goal is to estimate x from y.

The filter passes the frequencies where the signal dominates and attenuates those where the noise dominates. In practice f_y is estimated by Welch’s method, f_\varepsilon = \sigma^2 / (2\pi) for a white noise, and f_x = \max(f_y - f_\varepsilon, 0). The Wiener filter is optimal among linear filters only: it blurs the jumps of a piecewise-constant signal, whose energy spreads over all frequencies.

With the universal threshold \lambda = \sigma\sqrt{2\log n}, every pure-noise coefficient is set to zero with high probability, so the estimate does not invent oscillations. Donoho and Johnstone (1994) proved that its risk is within a factor of order 2\log n of an oracle that knows which coefficients to keep. The noise level is estimated from the finest wavelet details, \hat\sigma = \mathrm{median}(|d^{(1)}|) / 0.6745, which the few signal coefficients barely affect.

Soft thresholding shrinks every kept coefficient by \lambda, which biases large ones; hard thresholding keeps them intact but is discontinuous. The code compares both with a moving average on a piecewise-constant signal with a bump.

import numpy as np

def haar(x):                                 # orthonormal Haar transform (Prop. 1.12)
    coeffs, a = [], x.copy()
    while len(a) > 1:
        coeffs.append((a[0::2] - a[1::2]) / np.sqrt(2))   # details
        a = (a[0::2] + a[1::2]) / np.sqrt(2)              # approximation
    return a, coeffs

def ihaar(a, coeffs):
    for d in reversed(coeffs):
        out = np.empty(2 * len(a))
        out[0::2], out[1::2] = (a + d) / np.sqrt(2), (a - d) / np.sqrt(2)
        a = out
    return a

rng = np.random.default_rng(0)
n = 1024
t = np.arange(n) / n
clean = 2 * (t > 0.2) - 3 * (t > 0.45) + 1.5 * (t > 0.7) + np.exp(-((t - 0.85) / 0.03) ** 2)
y = clean + 0.4 * rng.standard_normal(n)

a, c = haar(y)
sigma = np.median(np.abs(c[0])) / 0.6745     # noise level from the finest details
lam = sigma * np.sqrt(2 * np.log(n))         # universal threshold (Prop. 6.4)
soft = lambda d: np.sign(d) * np.maximum(np.abs(d) - lam, 0)   # Theorem 6.3
hard = lambda d: d * (np.abs(d) > lam)                          # Proposition 6.5
x_soft = ihaar(a, [soft(d) for d in c])
x_hard = ihaar(a, [hard(d) for d in c])
x_ma = np.convolve(y, np.ones(9) / 9, mode="same")       # moving average, width 9
rmse = lambda u: np.sqrt(np.mean((u - clean) ** 2))
print(f"sigma estimate {sigma:.3f} (true 0.4), threshold {lam:.2f}")
for name, u in [("noisy", y), ("moving average", x_ma), ("Haar, soft", x_soft), ("Haar, hard", x_hard)]:
    print(f"RMSE {name:15s} {rmse(u):.3f}")

The estimated noise level is 0.376 for a true 0.4. The errors (RMSE) are 0.389 for the noisy signal, 0.167 for the moving average, 0.234 for soft thresholding and 0.131 for hard thresholding. Soft thresholding pays for its bias on the large jump coefficients; hard thresholding keeps the jumps sharp and wins.

For signals that are exactly piecewise constant, total-variation denoising solves \min_x \frac12\|y - x\|^2 + \lambda \sum_t |x_{t+1} - x_t| (Rudin, Osher and Fatemi, 1992). Its solution is piecewise constant, so it is a convex relative of the change-point detection of ch. 5; in one dimension it is computed exactly by a direct algorithm (Condat, 2013).

Missing values

Splines and linear interpolation are smooth and cheap, but they ignore the dynamics: across a long gap they draw a smooth curve, not a plausible piece of the series. A model of the dynamics does better.

Under MAR, m_{t \mid n} is the best imputation in mean square for a linear Gaussian model, and P_{t \mid n} says how uncertain it is. The code imputes an AR(2) series with scattered missing values and one gap of 20 samples.

import numpy as np

rng = np.random.default_rng(9)
n, phi = 1000, np.array([1.6, -0.8])         # quasi-periodic AR(2), period about 13
x = np.zeros(n)
for t in range(2, n):
    x[t] = phi @ x[t - 2 : t][::-1] + rng.standard_normal()
miss = np.zeros(n, bool)
miss[700:720] = True                         # a gap of 20 samples
miss[500:] |= rng.random(n - 500) < 0.15     # + 15% missing at random
y = np.where(miss, np.nan, x)

def kalman_smoother(y, F, Q, H, R):          # filter (Theorem 2.8) + RTS smoother
    n, d = len(y), F.shape[0]
    m, P = np.zeros(d), np.eye(d) * 10.0
    ms, Ps, mp, Pp = np.zeros((n, d)), np.zeros((n, d, d)), np.zeros((n, d)), np.zeros((n, d, d))
    for t in range(n):
        m, P = F @ m, F @ P @ F.T + Q        # predict
        mp[t], Pp[t] = m, P
        if not np.isnan(y[t]):               # update only when observed (Prop. 6.8)
            S = H @ P @ H.T + R
            G = P @ H.T / S
            m, P = m + (G * (y[t] - H @ m)).ravel(), P - G @ H @ P
        ms[t], Ps[t] = m, P
    for t in range(n - 2, -1, -1):           # backward pass
        J = Ps[t] @ F.T @ np.linalg.inv(Pp[t + 1])
        ms[t] = ms[t] + J @ (ms[t + 1] - mp[t + 1])
        Ps[t] = Ps[t] + J @ (Ps[t + 1] - Pp[t + 1]) @ J.T
    return ms[:, 0]

F = np.array([phi, [1.0, 0.0]])              # state (x_t, x_{t-1})
Q = np.diag([1.0, 0.0])
H = np.array([[1.0, 0.0]])
x_kal = kalman_smoother(y, F, Q, H, R=np.array([[1e-6]]))
obs = np.flatnonzero(~miss)
x_lin = np.interp(np.arange(n), obs, x[obs])
for name, u in [("linear interpolation", x_lin), ("AR(2) Kalman smoother", x_kal)]:
    gap, scattered = np.arange(700, 720), np.flatnonzero(miss & (np.arange(n) < 700))
    print(f"{name:22s} RMSE gap {np.sqrt(np.mean((u - x)[gap] ** 2)):.2f}"
          f" | scattered {np.sqrt(np.mean((u - x)[scattered] ** 2)):.2f}")
print("standard deviation of the series:", round(x.std(), 2))

On the scattered missing values, the smoother’s error is 0.46 against 0.79 for linear interpolation. In the 20-sample gap both errors (3.50 and 3.63) are close to the standard deviation of the series, 3.61. The gap is longer than the memory of the process, so no method can know the values in its middle; the smoother at least reports this through its variance P_{t \mid n}.

When several channels are recorded, they often share a few latent signals. The matrix of channels \times time is then nearly low rank, and the observed channels inform the missing ones.

import numpy as np

rng = np.random.default_rng(10)
T, d, r = 400, 20, 3                         # 20 sensors driven by 3 latent signals
t = np.arange(T)
latent = np.stack([np.sin(2 * np.pi * t / 50), np.sign(np.sin(2 * np.pi * t / 120)),
                   np.cumsum(rng.standard_normal(T)) / 10])
X = rng.standard_normal((d, r)) @ latent + 0.1 * rng.standard_normal((d, T))
mask = rng.random((d, T)) > 0.3              # 30% of the entries are missing

def soft_impute(Y, mask, lam, iters=200):    # Mazumder, Hastie and Tibshirani (2010)
    Z = np.zeros_like(Y)
    for _ in range(iters):
        U, s, Vt = np.linalg.svd(np.where(mask, Y, Z), full_matrices=False)
        Z = (U * np.maximum(s - lam, 0)) @ Vt    # singular-value soft thresholding
    return Z

Z = soft_impute(X, mask, lam=2.0)
lin = np.array([np.interp(t, t[m], row[m]) for row, m in zip(X, mask)])
err = lambda A: np.sqrt(np.mean((A - X)[~mask] ** 2))
print(f"RMSE on missing entries: per-sensor linear interpolation {err(lin):.3f}, "
      f"low-rank soft-impute {err(Z):.3f}")
print("top singular values:", [round(float(v), 1) for v in np.linalg.svd(Z, compute_uv=False)[:5]])

Soft-impute reaches an error of 0.140 against 0.256 for per-sensor interpolation, and the completed matrix has three dominant singular values (440.4, 92.6, 51.5) followed by values near 0: the rank 3 of the generating model is recovered. Chapter 8 adds a graph between the sensors to the same idea.

7. Extracting and selecting relevant features

Many learning problems on series are solved by turning each series into a vector of well-chosen numbers — its features — and then using any standard classifier or regressor. This chapter lists the main families of features, presents random convolutional features, and then selects the relevant features without fooling oneself. Named features also serve the official emphasis on understanding and interpreting the data.

Feature extraction

Large catalogues exist. The hctsa project of Fulcher and Jones (2014) computes thousands of features. catch22 (Lubba and coauthors, 2019) keeps 22 of them, chosen for accuracy and low redundancy across many data sets. tsfresh (Christ and coauthors, 2018) filters a large catalogue by hypothesis tests.

The dilations cover all scales, and the proportion of positive values measures how often a pattern is present, wherever it is. With K = 10\,000 kernels and a ridge classifier, ROCKET matched the most accurate classifiers of its time on the UCR archive at a small fraction of their cost. The block below is a simplified version (no padding, 500 kernels).

import numpy as np
from sklearn.linear_model import RidgeClassifierCV

rng = np.random.default_rng(12)
def cbf(n, length=128):                      # Cylinder-Bell-Funnel, as in ch. 3
    t, X, y = np.arange(length), [], np.repeat([0, 1, 2], n)
    for label in y:
        a = rng.integers(16, 33); b = a + rng.integers(32, 97)
        chi = ((t >= a) & (t <= b)).astype(float)
        shape = [chi, chi * (t - a) / (b - a), chi * (b - t) / (b - a)][label]
        X.append((6 + rng.standard_normal()) * shape + rng.standard_normal(length))
    return np.array(X), y

def random_kernels(k, T):                    # Definition 7.2
    ks = []
    for _ in range(k):
        L = rng.choice([7, 9, 11])
        w = rng.standard_normal(L); w -= w.mean()
        dil = int(2 ** rng.uniform(0, np.log2((T - 1) / (L - 1))))
        ks.append((w, rng.uniform(-1, 1), dil))
    return ks

def transform(X, ks):                        # 2 features per kernel: PPV and max
    out = []
    for w, b, dil in ks:
        wd = np.zeros((len(w) - 1) * dil + 1); wd[::dil] = w     # dilated kernel
        C = np.array([np.convolve(x, wd[::-1], mode="valid") for x in X]) + b
        out += [(C > 0).mean(1), C.max(1)]
    return np.column_stack(out)

Xtr, ytr = cbf(10); Xte, yte = cbf(100)
ks = random_kernels(500, Xtr.shape[1])
clf = RidgeClassifierCV(alphas=np.logspace(-3, 3, 10)).fit(transform(Xtr, ks), ytr)
print("random-kernel features: test accuracy", round(clf.score(transform(Xte, ks), yte), 3))

From 30 training series, the random features reach a test accuracy of 0.97 on 300 series, without any distance computation between series.

Feature selection

Extracting hundreds of features is easy; many are irrelevant or redundant. Selecting a few improves accuracy when data are scarce, speeds up the model, and makes it interpretable. The official survey of Li and coauthors (2018) organises the methods.

Mutual information thus extends the squared correlation, and it also detects nonlinear dependence. It is estimated by nearest-neighbour methods (Kraskov, Stögbauer and Grassberger, 2004), which scikit-learn uses. Ranking by I(X_j; Y) alone keeps redundant features: the mRMR criterion (Peng, Long and Ding, 2005) adds features greedily by relevance to Y minus mean redundancy with the features already chosen.

A feature enters the model only when its correlation with the residual reaches \lambda. For \lambda \ge \|X^\top y\|_\infty / n the solution is \hat\beta = 0; decreasing \lambda lets features enter one by one along the regularisation path, which ranks them. Among strongly correlated features the lasso picks one somewhat arbitrarily; the elastic net and stability selection (Meinshausen and Bühlmann, 2010) make the choice more stable.

The code builds two classes of AR(1) series, one of which contains a weak burst near 0.22 cycle per sample. It extracts 20 named features, adds 30 pure-noise features, and selects with mutual information inside the cross-validation and with an \ell_1-penalised logistic regression.

import numpy as np
from scipy.signal import lfilter, welch
from scipy.stats import skew, kurtosis
from sklearn.feature_selection import SelectKBest, mutual_info_classif
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

rng = np.random.default_rng(11)
N, T = 300, 256
y = rng.integers(0, 2, N)
X = np.zeros((N, T))
for k in range(N):                           # AR(1) background, random coefficient
    X[k] = lfilter([1.0], [1.0, -rng.uniform(0.2, 0.8)], rng.standard_normal(T))
    if y[k]:                                 # class 1: a weak burst near 0.22 cycle/sample
        f0, start = rng.uniform(0.2, 0.25), rng.integers(0, T - 64)
        X[k, start : start + 64] += 1.2 * np.sin(2 * np.pi * f0 * np.arange(64))

def features(x):                             # 20 named, interpretable features
    xc = x - x.mean()
    acf = [xc[:-h] @ xc[h:] / (xc @ xc) for h in range(1, 6)]
    f, P = welch(x, nperseg=64)
    bands = [P[(f >= lo) & (f < lo + 0.05)].sum() / P.sum() for lo in np.arange(0, 0.5, 0.05)]
    zc = np.mean(np.diff(np.sign(xc)) != 0)
    return np.r_[x.std(), skew(x), kurtosis(x), zc, np.max(np.abs(x)), acf, bands]

names = ["std", "skew", "kurt", "zero-cross", "max|x|"] + [f"acf{h}" for h in range(1, 6)] \
        + [f"band{lo:.2f}" for lo in np.arange(0, 0.5, 0.05)]
Phi = np.array([features(x) for x in X])
Phi = np.c_[Phi, rng.standard_normal((N, 30))]          # + 30 pure-noise features
names += [f"noise{j}" for j in range(30)]

mi = mutual_info_classif(Phi, y, random_state=0)
print("top-5 by mutual information:", [names[j] for j in np.argsort(mi)[::-1][:5]])
for k in (5, 50):                            # selection INSIDE the cross-validation
    pipe = make_pipeline(StandardScaler(), SelectKBest(mutual_info_classif, k=k),
                         LogisticRegression(max_iter=1000))
    print(f"k = {k:2d} features, 5-fold accuracy {cross_val_score(pipe, Phi, y, cv=5).mean():.3f}")
l1 = make_pipeline(StandardScaler(), LogisticRegression(l1_ratio=1.0, C=0.1, solver="liblinear"))
coef = l1.fit(Phi, y)[-1].coef_.ravel()
print("L1-logistic keeps:", [names[j] for j in np.flatnonzero(coef)])

Mutual information ranks first the power in the band 0.20–0.25, where the burst lives, then the autocorrelations at lags 2 to 4, which the oscillation also changes. One noise feature reaches the top five: mutual-information estimates are noisy on 300 samples. Five selected features give a cross-validated accuracy of 0.953, against 0.900 with all 50.

The \ell_1 model keeps four features: the standard deviation, the autocorrelation at lag 4 and the two bands 0.20 and 0.25. Each is a statement a domain expert can check: “class 1 has more energy between 0.20 and 0.30 cycle per sample”.

8. Multivariate series and graphs

With d sensors recorded together — electrodes on a scalp, stations of a weather network, joints of a skeleton — the series is a d \times T array. Treating each channel separately wastes the structure between them; a full multivariate model has too many parameters. A graph between the channels is the compromise: it says which channels are expected to behave alike. This chapter follows the graph signal processing of Ortega and coauthors (2018), in the official references.

From vector autoregressions to graphs

A VAR(p) has p d^2 parameters, too many when d is large: estimating it with a lasso penalty on each row of A learns a sparse directed graph. Most of this chapter uses instead an undirected graph that encodes similarity between channels.

Signals on a graph

The identity follows by expanding: x^\top D x = \sum_i d_i x_i^2 and x^\top W x = \sum_{i,j} W_{ij} x_i x_j, and each pair contributes W_{ij}(x_i^2 + x_j^2 - 2x_ix_j)/2 twice.

A signal is thus smooth exactly when its energy sits in the low graph frequencies. The next definition filters signals by frequency.

In Bayesian terms, \hat x is the posterior mean when x has the Gaussian prior with precision proportional to L (made invertible by a small multiple of the identity) and the observation noise is white. Interpolating a dead sensor from its neighbours is the graph analogue of the splines of ch. 6.

Learning the graph

Often no graph is given, and it must be learned from the multivariate series itself.

The graphical lasso (Friedman, Hastie and Tibshirani, 2008) estimates a sparse precision matrix from the empirical covariance S: \max_{\Theta \succ 0} \ \log\det\Theta - \mathrm{tr}(S\Theta) - \alpha \sum_{i \ne j} |\Theta_{ij}| .

If the snapshots are smooth signals with precision L + \varepsilon I — a Gaussian Markov random field on the graph — the support of the precision is exactly the set of edges.

Learning a graph on which the data are smooth is therefore a problem that is linear in the weights. Kalofolias (2016) solves \min_{W \ge 0,\ W = W^\top,\ \mathrm{diag}(W) = 0} \ \sum_{i,j} W_{ij} Z_{ij} - \alpha\, \mathbf{1}^\top \log(W\mathbf{1}) + \beta \|W\|_F^2 , where the logarithmic barrier prevents isolated nodes and the Frobenius term controls sparsity. Dong and coauthors (2016) learn the Laplacian directly under a trace constraint.

The block below places 30 sensors at random in the unit square, connects those closer than 0.3, and draws 2000 smooth snapshots with precision L + 0.1 I. It then interpolates 30 % of missing readings with Theorem 8.6, and learns the graph back with the graphical lasso.

import numpy as np
import networkx as nx
from sklearn.covariance import GraphicalLasso

rng = np.random.default_rng(13)
G = nx.random_geometric_graph(30, 0.3, seed=13)          # 30 sensors in the unit square
Lap = nx.laplacian_matrix(G).toarray().astype(float)
lam, U = np.linalg.eigh(Lap)                             # graph Fourier basis (Def. 8.3)

T = 2000                                                 # smooth signals: precision L + 0.1 I
X = U @ (rng.standard_normal((30, T)) / np.sqrt(lam + 0.1)[:, None])
Y = X + 0.3 * rng.standard_normal(X.shape)
print("share of energy in the 5 lowest graph frequencies:",
      round(float(((U.T @ X)[:5] ** 2).sum() / (X ** 2).sum()), 2))

mask = rng.random(Y.shape) > 0.3                         # 30% of readings missing
gamma = 0.1                                              # Tikhonov (Theorem 8.6)
Xhat = np.column_stack([np.linalg.solve(np.diag(m * 1.0) + gamma * Lap, m * y)
                        for m, y in zip(mask.T, Y.T)])
mean_imp = np.where(mask, Y, (Y * mask).sum(0) / mask.sum(0))   # fill with snapshot mean
err = lambda A: np.sqrt(np.mean((A - X)[~mask] ** 2))
print(f"RMSE on missing readings: snapshot mean {err(mean_imp):.3f}, graph Tikhonov {err(Xhat):.3f}")

prec = GraphicalLasso(alpha=0.02, max_iter=200).fit(X.T).precision_   # learn the graph back
iu = np.triu_indices(30, 1)
top = np.argsort(-np.abs(prec[iu]))[: G.number_of_edges()]
true = {(min(e), max(e)) for e in G.edges()}
hits = sum((int(iu[0][k]), int(iu[1][k])) in true for k in top)
print(f"graphical lasso: {hits} of its top {G.number_of_edges()} entries are true edges")

The five lowest graph frequencies carry 83 % of the energy: the snapshots are smooth. On the missing readings, the graph interpolation reaches an error of 0.416 against 0.570 for the mean of the observed sensors. The graphical lasso, given 2000 snapshots, ranks the 114 true edges exactly first: 114 \text{ of the top } 114 \text{ entries are edges}.

With fewer snapshots, or with signals that are only approximately smooth, the recovery degrades, and the choice of \alpha becomes the main difficulty.

9. Neural networks for time series

Self-study chapter, beyond the official syllabus. Neural networks learn the representation instead of choosing it. This chapter gives the three architectures used on series, one fact about each, and the baselines they must be compared with; Deep Learning covers the general theory.

With a nonlinearity, each factor W is multiplied by a diagonal matrix of derivatives, which does not change the picture. The LSTM (Hochreiter and Schmidhuber, 1997) adds a gated memory cell whose self-connection is close to the identity, so that information and gradients can travel over long lags. DeepAR (Salinas and coauthors, 2020) trains one recurrent network on many related series and outputs a predictive distribution.

The block trains a three-layer dilated CNN (kernel 7, dilations 1, 2, 4) on 300 Cylinder–Bell–Funnel series, with global max pooling so that the class does not depend on the position of the pattern.

Its convolutions are centred, with padding on both sides, rather than causal as in Definition 9.1. For classification the whole series is available, so causality is not needed; Proposition 9.2 still gives the receptive field, now centred on t.

import numpy as np
import torch
from torch import nn

rng = np.random.default_rng(14)
torch.manual_seed(0)
torch.set_num_threads(1)                     # small model: one thread is fastest
def cbf(n, length=128):                      # Cylinder-Bell-Funnel, as in ch. 3
    t, X, y = np.arange(length), [], np.repeat([0, 1, 2], n)
    for label in y:
        a = rng.integers(16, 33); b = a + rng.integers(32, 97)
        chi = ((t >= a) & (t <= b)).astype(float)
        shape = [chi, chi * (t - a) / (b - a), chi * (b - t) / (b - a)][label]
        X.append((6 + rng.standard_normal()) * shape + rng.standard_normal(length))
    return torch.tensor(np.array(X), dtype=torch.float32)[:, None, :], torch.tensor(y)

Xtr, ytr = cbf(100)
Xte, yte = cbf(100)
net = nn.Sequential(                         # dilated convolutions, then global pooling
    nn.Conv1d(1, 16, 7, padding=3), nn.ReLU(),
    nn.Conv1d(16, 16, 7, padding=6, dilation=2), nn.ReLU(),
    nn.Conv1d(16, 16, 7, padding=12, dilation=4), nn.ReLU(),
    nn.AdaptiveMaxPool1d(1), nn.Flatten(), nn.Linear(16, 3))
opt = torch.optim.Adam(net.parameters(), lr=3e-3)
for epoch in range(20):
    for idx in torch.randperm(len(ytr)).split(32):
        opt.zero_grad()
        loss = nn.functional.cross_entropy(net(Xtr[idx]), ytr[idx])
        loss.backward()
        opt.step()
with torch.no_grad():
    acc = (net(Xte).argmax(1) == yte).float().mean().item()
print(f"1-D CNN, receptive field {1 + 6 * (1 + 2 + 4)} samples: test accuracy {acc:.3f}")

The receptive field is 1 + 6 \times 7 = 43 samples (Proposition 9.2), and the network classifies all 300 test series correctly. It needed 300 training series; 1-NN DTW and the random kernels of ch. 7 reached similar accuracy from 30. Exercise 15 compares the three methods as the training set shrinks.

10. Research directions and open problems

Self-study chapter. Each theme of the course is an active research field. This chapter points to a few directions, with works whose authors and dates are given so that they can be found; it is a starting point for the mini-project, not a survey.

Foundation models for forecasting

Following language models, several groups pretrain one large model on many series and use it without retraining on a new one (“zero-shot”). Chronos (Ansari and coauthors, 2024) scales each series, quantises its values into a vocabulary of bins, and trains a language-model architecture on the resulting tokens. TimesFM (Das, Kong, Sen and Zhou, 2024) is a decoder-only model pretrained on a large corpus of real and synthetic series.

Open questions: when do such models beat a well-tuned local model fitted to the series at hand? How to evaluate them when public benchmarks may be part of the pretraining corpus? Are their predictive intervals calibrated?

Honest evaluation

Benchmarks drive the field, and several have been questioned. For anomaly detection, Wu and Keogh (2023) found trivial and mislabelled anomalies in popular data sets, and Kim and coauthors (2022) showed that the point-adjust protocol inflates scores (Exercise 16). For classification, the UCR archive (Dau and coauthors, 2019) and the UEA multivariate archive (Bagnall and coauthors, 2018) standardise the splits, and comparisons over many data sets call for rank-based statistical tests (Demšar, 2006).

Change points beyond the mean

The official review of Truong, Oudre and Vayatis (2020) organises the field by cost, search method and constraint. Kernel costs detect changes in distribution, but the choice of the kernel and of the penalty remains largely empirical (Arlot, Celisse and Harchaoui, 2019). Online detection with false-alarm guarantees, change points in high-dimensional series, and costs learned from annotated examples are active directions.

Irregular sampling and missing values

Clinical and environmental series are sampled irregularly, with informative missingness. GRU-D (Che and coauthors, 2018) feeds the missingness mask and the time since the last observation into a recurrent network. Latent ordinary differential equations (Rubanova, Chen and Duvenaud, 2019) model the hidden state in continuous time. CSDI (Tashiro and coauthors, 2021) imputes with a conditional diffusion model and returns samples, hence uncertainty. Imputation under MNAR mechanisms, with calibrated uncertainty, is largely open.

Interpretable pattern extraction

Convolutional dictionary learning gives atoms that experts can inspect. On electrophysiological signals, rank-one multivariate atoms (Dupré la Tour and coauthors, 2018) and heavy-tailed noise models (Jas and coauthors, 2017) produce waveforms that match known physiology. Scaling these methods to days of recordings, and choosing the number and length of the atoms automatically, are open problems. Shapelet and random-kernel classifiers raise the same question: how to explain a decision with a few patterns without losing accuracy?

Graphs that change in time

Spatio-temporal neural networks combine graph filters with temporal models: DCRNN (Li and coauthors, 2018) uses diffusion convolutions on a road network inside a recurrent network, and Graph WaveNet (Wu and coauthors, 2019) learns an adaptive adjacency matrix together with dilated temporal convolutions. Time–vertex signal processing (Grassi and coauthors, 2018) gives the linear theory. Learning graphs that change over time from few samples, with guarantees, is open.

Open problems to think about

• A principled penalty for change-point detection under dependent, non-Gaussian noise.
• Evaluation protocols for anomaly detection that do not depend on an arbitrary threshold or event definition.
• Imputation under MNAR mechanisms with calibrated uncertainty.
• Zero-shot forecasting: conditions under which pretraining on other series provably helps.
• Interpretable models that match deep networks on classification accuracy.

Exercises

Levels: ★ apply the course, ★★ combine results or prove a short statement, ★★★ open-ended or research-flavoured. The Practice column of the session plan says when to do each one.

1. ★ The MA(1) autocorrelation. Let X_t = \varepsilon_t + \theta\varepsilon_{t-1}, where \varepsilon is a white noise of variance \sigma^2. Compute the autocovariance at every lag, show that |\rho(1)| \le 1/2, and find the values of \theta for which \rho(1) = 0.4.

   Solution: \gamma(0) = \sigma^2(1 + \theta^2), \gamma(\pm 1) = \theta\sigma^2 and \gamma(h) = 0 for |h| \ge 2, so \rho(1) = \theta / (1 + \theta^2).

   Since (1 - |\theta|)^2 \ge 0, 1 + \theta^2 \ge 2|\theta| and |\rho(1)| \le 1/2, with equality at \theta = \pm 1.

   The equation 0.4\theta^2 - \theta + 0.4 = 0 gives \theta = 2 or \theta = 0.5: two models with the same autocorrelation. One keeps the invertible one, |\theta| < 1, whose shocks can be recovered from the past of X.

2. ★ The random walk is not stationary. Let S_t = \varepsilon_1 + \dots + \varepsilon_t. Compute \mathrm{Cov}(S_t, S_{t+h}) and \mathrm{Corr}(S_t, S_{t+h}) for h \ge 0. Explain why the sample ACF of a random walk decays very slowly, and what differencing does.

   Solution: \mathrm{Cov}(S_t, S_{t+h}) = t\sigma^2 (the common shocks), so \mathrm{Corr}(S_t, S_{t+h}) = \sqrt{t / (t + h)}, which depends on t (no stationarity) and is close to 1 for h \ll t. The sample ACF therefore stays near 1 for many lags. The differenced series \nabla S_t = \varepsilon_t is a white noise (Example 1.6).

3. ★ Finding a hidden period (Python). Generate x_t = \cos(2\pi \cdot 0.13\, t) + 2\varepsilon_t for t < 1000 with Gaussian \varepsilon. Locate the peak of the raw periodogram and of Welch’s estimate, and compare their fluctuations at frequencies above 0.3, where there is only noise.

   Hints and expected result: both peaks are at 0.13. Off the peak, the raw periodogram has a relative spread (standard deviation over mean) of about 1, as predicted by the \chi^2_2 law of Proposition 1.9; Welch’s estimate has about 0.33.

   import numpy as np
   from scipy.signal import periodogram, welch
   
   rng = np.random.default_rng(0)
   t = np.arange(1000)
   x = np.cos(2 * np.pi * 0.13 * t) + 2 * rng.standard_normal(1000)
   f1, P1 = periodogram(x)                      # raw periodogram
   f2, P2 = welch(x, nperseg=200)               # Welch, segments of 200 samples
   print("peak frequencies:", f1[P1.argmax()], f2[P2.argmax()])
   far1, far2 = P1[f1 > 0.3], P2[f2 > 0.3]      # noise-only frequencies
   print("relative spread off the peak: raw", round(far1.std() / far1.mean(), 2),
         "| Welch", round(far2.std() / far2.mean(), 2))

4. ★ Forecast intervals of an AR(1). For \phi = 0.6, \sigma = 0.5 and the last value x_t = 2, give the forecasts at horizons 1, 2 and 3 and their Gaussian 95 % intervals. What is the limit of the interval width?

   Solution: by Example 2.5, the forecasts are 1.2, 0.72 and 0.432. The error variances are 0.25, 0.34 and 0.3724, that is, \sigma^2 times 1, 1.36 and 1.4896. The standard deviations 0.5, 0.583 and 0.610, multiplied by 1.96, give the half-widths 0.98, 1.14 and 1.20.

   The intervals are [0.22, 2.18], [-0.42, 1.86] and [-0.76, 1.63]. The variance tends to \sigma^2 / (1 - \phi^2) = 0.39, so the half-width tends to 1.96 \times 0.625 = 1.23.

5. ★ DTW by hand. Compute \mathrm{DTW}(x, y) for x = (1, 3, 4) and y = (1, 2, 3, 4) with the recursion of Theorem 3.3, and give an optimal path.

   Solution: the table D(i, j), rows x_i and columns y_j, is \begin{pmatrix} 0 & 1 & 5 & 14 \\ 4 & 1 & 1 & 2 \\ 13 & 5 & 2 & 1 \end{pmatrix}, so \mathrm{DTW}(x, y)^2 = 1 and \mathrm{DTW}(x, y) = 1.

   The path (1,1), (2,2), (2,3), (3,4) pays only for matching 3 with 2. The path (1,1), (1,2), (2,3), (3,4), which matches 1 with 2, is also optimal: optimal paths need not be unique.

6. ★ Seasonal differencing. A monthly series is x_t = 10 + 0.2\,t + s_t + \varepsilon_t, with s of period 12 and \varepsilon a white noise. What is \nabla_{12} x_t? What is its autocorrelation at lag 12? What does one more ordinary difference do?

   Solution: by Proposition 1.13, \nabla_{12} x_t = 2.4 + \varepsilon_t - \varepsilon_{t-12}, a constant mean 12 \times 0.2 plus an MA(12) noise with coefficient -1 at lag 12. Its autocorrelation at lag 12 is -1/2, and 0 at the other nonzero lags. One more difference removes the constant 2.4 and leaves a zero-mean stationary series. Differencing twice at lag 1 instead would not remove the season.

7. ★★ Yule–Walker for an AR(2). A stationary series has \gamma(0) = 1, \rho(1) = 0.5 and \rho(2) = 0.1. Fit an AR(2) by Theorem 2.3, give \sigma^2, check that the model is causal, and predict \rho(3).

   Solution: the equations 0.5 = \phi_1 + 0.5\phi_2 and 0.1 = 0.5\phi_1 + \phi_2 give \phi_2 = -0.2 and \phi_1 = 0.6. Then \sigma^2 = 1 - 0.6 \times 0.5 + 0.2 \times 0.1 = 0.72.

   The roots of 1 - 0.6z + 0.2z^2 are complex with |z|^2 = 1/0.2 = 5, so |z| \approx 2.24 > 1: the model is causal. Finally \rho(3) = 0.6 \times 0.1 - 0.2 \times 0.5 = -0.04.

8. ★★ Soft thresholding and the fixed point of ISTA. (a) Compute S_1(u) for u = (3, -0.5, 1.2, -2). (b) Show that z^* minimises the lasso F of Proposition 4.5 if and only if it is a fixed point of the ISTA iteration.

   Solution: (a) S_1(u) = (2, 0, 0.2, -1).

   2. z^* is optimal iff 0 \in D^\top(Dz^* - y) + \lambda\,\partial\|z^*\|_1. Multiply by 1/L and set v = z^* - \frac1L D^\top(Dz^* - y): the condition reads v - z^* \in \frac{\lambda}{L}\partial\|z^*\|_1.

   This is exactly the optimality condition of the proximal problem whose solution is S_{\lambda/L}(v). Hence z^* = S_{\lambda/L}(v), the ISTA update applied to z^*.

9. ★★ Correlation, z-normalisation and a lower bound. (a) Deduce from Proposition 3.10 that 1-NN with the z-normalised Euclidean distance is 1-NN with the largest Pearson correlation. (b) For x = (0, 0, 3, 0), y = (0, 1, 1, 0) and w = 1, compute the envelope of y, \mathrm{LB}(x, y) and \mathrm{DTW}_w(x, y).

   Solution: (a) d^2 = 2m(1 - \rho) is a decreasing function of \rho, so both rules select the same neighbour.

   2. U = (1, 1, 1, 1) and L = (0, 0, 0, 0). Only x_3 = 3 leaves the envelope, so \mathrm{LB}^2 = (3 - 1)^2 = 4 and \mathrm{LB} = 2.

   For the banded DTW, the value 3 of x must be matched with one of y_2, y_3, y_4, at cost at least 4. Every sample of y is matched too, so the other value 1 of y is matched either with 3 (cost 4) or with a 0 of x (cost 1).

   The minimum is \mathrm{DTW}_w^2 = 4 + 1 = 5, and indeed \mathrm{LB} = 2 \le \sqrt 5 \approx 2.24.

10. ★★ Spikes and sinusoids. Let D = [I, F] in dimension n, where F is the orthonormal DFT matrix with entries e^{-2i\pi kt/n}/\sqrt n.

    1. Show that \mu(D) = 1/\sqrt n. (b) What sparsity does Theorem 4.4 guarantee for n = 64? (c) For n = m^2, show that the Dirac comb with m spikes spaced by m has two different representations with m atoms each, and conclude about the order of the bound.

    Solution: (a) Columns within I, or within F, are orthogonal; |\langle e_t, f_k \rangle| = 1/\sqrt n for every pair.

    2. \mu = 1/8, so uniqueness and recovery hold for \|x\|_0 < 4.5, that is, up to 4 atoms.

    3. The DFT of the comb \sum_{j < m} e_{jm} is again a comb of m spikes, so the signal is a sum of m spikes and also of m Fourier atoms.

    Hence \mathrm{spark}(D) \le 2m = 2\sqrt n, and uniqueness fails at sparsity \sqrt n. The coherence bound (1 + \sqrt n)/2 is therefore sharp up to a factor 2.

11. ★★ Binary segmentation misses a short segment (Python). Build 300 samples of \mathcal{N}(0, 1), 20 samples of \mathcal{N}(1.5, 1) and 300 more of \mathcal{N}(0, 1). Run binary segmentation (Note 5.7) and optimal partitioning (Proposition 5.5) with c_{L_2} and \beta = 2\log n. Explain the difference with Proposition 5.3.

    Hints and expected result: optimal partitioning finds changes at 299 and 320; binary segmentation finds none.

    The best single split leaves the short bump inside a long segment on either side: its gain is about \frac{300 \times 320}{620}(0.09)^2 \approx 1.4, far below \beta \approx 12.9, so the recursion stops at once. The pair of splits reduces the cost by about 20 \times 1.5^2 \approx 45 > 2\beta, which the exact search sees.

    import numpy as np
    
    rng = np.random.default_rng(0)
    y = np.r_[rng.standard_normal(300), 1.5 + rng.standard_normal(20), rng.standard_normal(300)]
    n = len(y)
    S1, S2 = np.r_[0, np.cumsum(y)], np.r_[0, np.cumsum(y**2)]
    cost = lambda a, b: S2[b] - S2[a] - (S1[b] - S1[a]) ** 2 / (b - a)
    beta = 2 * np.log(n)
    
    def binseg(a, b, min_size=5):                # Note 5.7: best single split, then recurse
        gains = [(cost(a, b) - cost(a, t) - cost(t, b), t) for t in range(a + min_size, b - min_size)]
        if not gains:
            return []
        g, t = max(gains)
        return binseg(a, t) + [t] + binseg(t, b) if g > beta else []
    
    def opt_partition(min_size=5):               # Proposition 5.5 (exact)
        F, last = np.full(n + 1, np.inf), np.zeros(n + 1, int)
        F[0] = -beta
        for t in range(min_size, n + 1):
            vals = [F[s] + cost(s, t) + beta for s in range(0, t - min_size + 1)]
            last[t] = int(np.argmin(vals)); F[t] = vals[last[t]]
        cps, t = [], n
        while t > 0:
            t = last[t]; cps.append(int(t))
        return sorted(cps)[1:]
    
    print("binary segmentation:", binseg(0, n))
    print("optimal partitioning:", opt_partition())

12. ★★ Exponential smoothing is a Kalman filter. In the local level model of Definition 2.7, with noise variances q and r, show that the one-step forecast follows \hat y_{t+1} = \hat y_t + K_t(y_t - \hat y_t). Find the steady-state gain for q = 1 and r = 4, and say what happens when y_t is missing.

    Solution: here F = H = 1, so the forecast \hat y_{t+1} = m_{t+1 \mid t} equals m_{t \mid t}, and m_{t+1 \mid t} = m_{t \mid t-1} + K_t\,(y_t - m_{t \mid t-1}): simple exponential smoothing with \alpha = K_t (Note 2.6).

    The predicted variance P satisfies P_{t+1 \mid t} = P_{t \mid t-1} r / (P_{t \mid t-1} + r) + q. At the fixed point, P^2 - qP - qr = 0, so P = (q + \sqrt{q^2 + 4qr})/2.

    With q = 1 and r = 4: P = (1 + \sqrt{17})/2 \approx 2.56 and K = P/(P + r) \approx 0.39.

    When y_t is missing, K_t = 0 (Proposition 6.8): the forecast stays put and its variance grows by q.

13. ★★ Selection bias in numbers (Python). Draw 50 samples with 2000 pure-noise features and random binary labels. Estimate by 5-fold cross-validation the accuracy of a logistic regression on the 20 features with the best F-score, (a) selected once on all the data, (b) selected inside each training fold.

    Hints and expected result: (a) gives about 0.96, (b) about 0.58, close to the true 0.5 (Note 7.6).

    import numpy as np
    from sklearn.feature_selection import SelectKBest, f_classif
    from sklearn.linear_model import LogisticRegression
    from sklearn.model_selection import cross_val_score
    from sklearn.pipeline import make_pipeline
    
    rng = np.random.default_rng(0)
    X, y = rng.standard_normal((50, 2000)), rng.integers(0, 2, 50)   # pure noise
    sel = SelectKBest(f_classif, k=20).fit(X, y)                       # WRONG: on all data
    wrong = cross_val_score(LogisticRegression(), sel.transform(X), y, cv=5).mean()
    pipe = make_pipeline(SelectKBest(f_classif, k=20), LogisticRegression())
    right = cross_val_score(pipe, X, y, cv=5).mean()                   # selection in each fold
    print(f"selection outside CV: {wrong:.2f}   inside CV: {right:.2f}")

14. ★★ The path graph. Let G be the path on d nodes with unit weights. (a) For k = 0, \dots, d - 1, check that the vector u_k(i) = \cos\big(\pi k (i - \tfrac12)/d\big), \quad i = 1, \dots, d, is an eigenvector of L with eigenvalue \lambda_k = 2 - 2\cos(\pi k / d). (b) For d = 4 and \gamma = 1, give the response of the Tikhonov filter of Theorem 8.6 at each frequency.

    Solution: (a) At an interior node, (Lu)(i) = 2u(i) - u(i-1) - u(i+1). The identity \cos(a - \theta) + \cos(a + \theta) = 2\cos a\cos\theta, with \theta = \pi k/d, gives (2 - 2\cos\theta)\,u(i).

    At node 1, (Lu)(1) = u(1) - u(2). The formula extended to i = 0 gives u(0) = u(1), so the interior computation still applies; symmetrically at node d. These are the cosines of the DCT-II.

    2. The eigenvalues are 0, 2 - \sqrt 2 \approx 0.586, 2 and 2 + \sqrt 2 \approx 3.414. The responses 1/(1 + \lambda_k) are 1, 0.631, 0.333 and 0.227: a low-pass filter.

15. ★★★ Three classifiers and the size of the training set (Python). On Cylinder–Bell–Funnel, compare 1-NN DTW (ch. 3), the random-kernel features of ch. 7 and the CNN of ch. 9 with 5, 10, 30 and 100 training series per class, over 5 random draws. Add a harder variant in which each series is also randomly stretched in time by a factor in [0.8, 1.2]. Report accuracies, training and prediction times, and discuss.

    Hints and expected result: with 10 series per class, the code of ch. 3 gives no error for 1-NN DTW and that of ch. 7 an accuracy of 0.97 for the random kernels. Expect the CNN to need more data and more tuning, and to be the most sensitive to the random seed. 1-NN DTW is the slowest at prediction time (one DTW per training series); the random kernels are the fastest overall. Stretching barely affects DTW, which is built for it. A good report gives mean and spread over the draws, not a single number.

16. ★★★ How point-adjust inflates anomaly scores (Python, open-ended). On 10 000 samples with 5 anomalous segments of 100 samples, score every sample with a uniform random number. Compute the best F1 over thresholds point by point, then after point-adjust (a whole segment counts as detected if one of its points is flagged). Then propose an event-level metric that rewards early detection and penalises false alarms, and test it on a real detector of ch. 5.

    Hints and expected result: the best point-wise F1 of the random detector is about 0.09, and the best point-adjusted F1 is about 0.79: the inflation described by Kim and coauthors (2022). A sound metric counts events, not points, gives credit that decreases with the detection delay, and counts each false-alarm event once. Its value must stay near 0 for a random detector — this is the sanity check.

    import numpy as np
    
    rng = np.random.default_rng(0)
    n = 10_000
    label = np.zeros(n, bool)
    for s in rng.choice(np.arange(0, n - 100, 200), 5, replace=False):
        label[s : s + 100] = True                # 5 anomalous segments of 100 samples
    score = rng.random(n)                        # a detector that knows nothing
    
    def f1(pred, label):
        tp = np.sum(pred & label); fp = np.sum(pred & ~label); fn = np.sum(~pred & label)
        return 2 * tp / (2 * tp + fp + fn)
    
    def point_adjust(pred, label):               # a hit anywhere marks the whole segment
        pred = pred.copy()
        edges = np.flatnonzero(np.diff(np.r_[0, label.astype(int), 0]))
        for a, b in zip(edges[::2], edges[1::2]):
            if pred[a:b].any():
                pred[a:b] = True
        return pred
    
    ths = np.quantile(score, np.linspace(0.5, 0.9999, 200))
    print("best point-wise F1   :", round(max(f1(score > h, label) for h in ths), 3))
    print("best point-adjust F1 :", round(max(f1(point_adjust(score > h, label), label) for h in ths), 3))

Projects and assessment

The official assessment (master-mva.com) combines lab reports (25 %) and mini-projects graded on the report (25 %), the source code (25 %) and the oral defence (25 %); the Université Paris-Saclay course detail describes it as non-integral continuous assessment. No written exam is announced. The subjects below are sized as mini-projects with a report, code and a defence; the written exam is only a self-assessment proposal.

Solo projects

1. [Solo · ≈ 20 h] Paper study: the shapelet transform (Hills and coauthors, 2014). Question: when the class is defined by a short pattern at a random position, does the shapelet transform beat 1-NN DTW?
   • Implement the candidate search with an F-statistic quality, the removal of self-similar candidates, the transform and a linear classifier.
   • Compare with 1-NN DTW on Cylinder–Bell–Funnel and on the random-position data set of ch. 3, varying noise level and series length.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 3 and 7.
   • Success: test accuracies of both methods reported as mean and standard deviation over 10 random draws, for at least three training-set sizes, and the question answered by a paired test over the draws. As a sanity check, the shapelet transform must beat Euclidean 1-NN on the random-position data, as in ch. 3 (0.917 against 0.67).
   • Also required: the three best shapelets plotted next to the planted triangle and plateau patterns, with their correlation to them.
2. [Solo · ≈ 20 h] Paper study: PELT (Killick, Fearnhead and Eckley, 2012). Question: is PELT exact and linear in practice?
   • Implement optimal partitioning, PELT and binary segmentation for the quadratic, Gaussian and kernel costs (ch. 5).
   • Check that PELT and optimal partitioning return identical segmentations; plot the run time against n from 10^3 to 10^5 on a log–log scale.
   • Calibrate the penalty (BIC-type against an elbow rule) and report precision and recall with a margin.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Success: identical outputs and a slope close to 1 for PELT when the number of changes grows with n.
3. [Solo · ≈ 20 h] Linear or nonlinear denoising? Question: on which signals does wavelet thresholding beat the optimal linear filter?
   • Implement the Wiener filter with an estimated spectrum (Theorem 6.1), Haar soft and hard thresholding with the universal threshold, and their translation-invariant version by cycle spinning (Coifman and Donoho, 1995).
   • Test on piecewise-constant, smooth and chirp signals at three noise levels.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 1 and 6.
   • Success: a table of errors whose ranking is explained by Theorem 6.1 and Propositions 6.4–6.5.

Group projects

1. [Group of 2–3 · ≈ 30 h each] Paper study: convolutional dictionary learning (Wohlberg, 2015; Dupré la Tour and coauthors, 2018).
   • Implement CSC by FISTA and by ADMM in the Fourier domain (Proposition 4.9), then alternating dictionary learning.
   • Apply it to synthetic multichannel signals built from rank-one atoms, then optionally to a public MEG recording such as the MNE-Python sample data set.
   • Split: one member writes and benchmarks the two CSC solvers, one the dictionary update and its initialisation, one the multichannel rank-one extension and the evaluation.
   • Deliverables: code, a 6–8 page report, a 15-minute group talk. Chapters 0 and 4. Success: atoms recovered with correlation above 0.9 and activation F1 above 0.9 at a signal-to-noise ratio of 10 dB.
2. [Group of 2–3 · ≈ 30 h each] An honest anomaly-detection benchmark.
   • Generate series with labelled point, contextual and collective anomalies on a seasonal AR background (ch. 5).
   • Implement four detectors: seasonal residuals, matrix-profile discords (ch. 3), isolation forest on windows and CUSUM.
   • Evaluate with point-wise, point-adjusted and event-level metrics, and reproduce the inflation reported by Kim and coauthors (2022).
   • Split: data generator and metrics; detectors; analysis and report. Deliverables: code, a 6–8 page report, a 15-minute talk. Success: a random detector scores near 0 with the event-level metric, and a table ranks the detectors by anomaly type.
3. [Group of 3 · ≈ 30 h each] Imputing a sensor network with a learned graph.
   • Start from synthetic Gaussian Markov random fields on random geometric graphs (ch. 8), then use a public multi-sensor data set, for instance the METR-LA traffic data released with DCRNN (Li and coauthors, 2018).
   • Learn the graph (graphical lasso, smoothness-based learning), then impute with graph Tikhonov, a time–vertex regulariser, low-rank soft-impute and per-channel Kalman smoothing (ch. 2, 6, 8).
   • Compare random and block missingness. Split: graph learning; imputation methods; missingness scenarios and evaluation.
   • Deliverables: code, a 6–8 page report, a 15-minute talk. Success: at least 20 % lower error than per-channel linear interpolation under block missingness.

Written exam (proposal)

Format: 2 hours, one handwritten A4 sheet allowed, no electronic device, four exercises.

1. Question (4 points). Define weak stationarity and the spectral density of a stationary process. State the relation between the spectral density and the autocovariance. Explain why the raw periodogram is not a consistent estimator of the spectral density, and how Welch’s method helps. Expected: Definition 1.3 and Theorem 1.7; the periodogram of a Gaussian white noise follows a scaled \chi^2_2 law whose variance does not decrease with n (Proposition 1.9); averaging K segments divides the variance by about K and coarsens the frequency resolution by the same factor.

2. Question (5 points). A stationary series has \gamma(0) = 2, \rho(1) = 0.6 and \rho(2) = 0.2. Fit an AR(2) by the Yule–Walker equations, give the noise variance and check causality. Forecast the next two values from x_t = 1 and x_{t-1} = 0.5, with their error variances.

   Expected: \phi_1 = 0.75, \phi_2 = -0.25 and \sigma^2 = 2(1 - 0.45 + 0.05) = 1.2. The roots of 1 - 0.75z + 0.25z^2 are complex with modulus 2, so the model is causal.

   The forecasts are 0.625 and 0.75 \times 0.625 - 0.25 = 0.219, with error variances 1.2 and 1.2(1 + 0.75^2) = 1.875.

3. Question (5 points). Let y = (1.0, 1.2, 0.8, 3.1, 2.9, 3.0, 3.2) with noise level \sigma = 0.2. Using Proposition 5.3, find the best single change point, decide whether to keep it with the penalty \beta = 2\sigma^2\log n, and decide whether a second change point should be added.

   Expected: the gains for \tau = 1, \dots, 6 are 1.60, 3.21, 7.20, 3.90, 2.41 and 1.23, so \tau = 3; with \beta \approx 0.156 it is kept.

   Inside the two segments, the best gains are 0.06 and 0.03, both below \beta. There is one change point, between the segment means 1.0 and 3.05.

4. Question (6 points). (a) Prove that for an orthonormal matrix W, the minimiser of \frac12\|y - x\|^2 + \lambda\|Wx\|_1 is W^\top S_\lambda(Wy), including the computation of the scalar proximal operator (3 points).

   2. Using \mathbb{P}(|Z| > u) \le \sqrt{2/\pi}\, e^{-u^2/2}/u for a standard Gaussian Z, prove that the largest absolute value of n i.i.d. \mathcal{N}(0, \sigma^2) variables exceeds \sigma\sqrt{2\log n} with probability at most 1/\sqrt{\pi\log n} (3 points).

   Expected: (a) the change of variables \theta = Wx and the case analysis of Proposition 4.5, as in Theorem 6.3; (b) the union bound and e^{-\log n} = 1/n, as in Proposition 6.4.

Grading rubric (proposal)

The weights are those of the official assessment; the descriptions are a proposal.

References

Official references (master-mva.com)

• Wohlberg, B. (2015). Efficient algorithms for convolutional sparse representations. IEEE Transactions on Image Processing, 25(1), 301–315 (journal issue of January 2016).
• Hills, J., Lines, J., Baranauskas, E., Mapp, J., & Bagnall, A. (2014). Classification of time series by shapelet transformation. Data Mining and Knowledge Discovery, 28(4), 851–881.
• Bisgaard, S., & Kulahci, M. (2011). Time Series Analysis and Forecasting by Example. John Wiley & Sons.
• Tošić, I., & Frossard, P. (2011). Dictionary learning. IEEE Signal Processing Magazine, 28(2), 27–38.
• Truong, C., Oudre, L., & Vayatis, N. (2020). Selective review of offline change point detection methods. Signal Processing, 167, 107299.
• Chandola, V., Banerjee, A., & Kumar, V. (2009). Anomaly detection: A survey. ACM Computing Surveys, 41(3), 1–58.
• Li, J., Cheng, K., Wang, S., Morstatter, F., Trevino, R. P., Tang, J., & Liu, H. (2018). Feature selection: A data perspective. ACM Computing Surveys, 50(6), 94.
• Ortega, A., Frossard, P., Kovačević, J., Moura, J. M. F., & Vandergheynst, P. (2018). Graph signal processing: Overview, challenges, and applications. Proceedings of the IEEE, 106(5), 808–828.

Textbooks

• Brockwell, P. J., & Davis, R. A. (1991). Time Series: Theory and Methods (2nd ed.). Springer.
• Box, G. E. P., Jenkins, G. M., Reinsel, G. C., & Ljung, G. M. (2015). Time Series Analysis: Forecasting and Control (5th ed.). Wiley.
• Hyndman, R. J., & Athanasopoulos, G. (2021). Forecasting: Principles and Practice (3rd ed.). OTexts.
• Durbin, J., & Koopman, S. J. (2012). Time Series Analysis by State Space Methods (2nd ed.). Oxford University Press.
• Mallat, S. (2009). A Wavelet Tour of Signal Processing: The Sparse Way (3rd ed.). Academic Press.
• Elad, M. (2010). Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing. Springer.
• Little, R. J. A., & Rubin, D. B. (2019). Statistical Analysis with Missing Data (3rd ed.). Wiley.

Representations, models and forecasting (ch. 1–2)

• Welch, P. D. (1967). The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms. IEEE Transactions on Audio and Electroacoustics, 15(2), 70–73.
• Cleveland, R. B., Cleveland, W. S., McRae, J. E., & Terpenning, I. (1990). STL: A seasonal-trend decomposition procedure based on loess. Journal of Official Statistics, 6(1), 3–73.
• Kalman, R. E. (1960). A new approach to linear filtering and prediction problems. Journal of Basic Engineering, 82(1), 35–45.
• Rauch, H. E., Tung, F., & Striebel, C. T. (1965). Maximum likelihood estimates of linear dynamic systems. AIAA Journal, 3(8), 1445–1450.
• Hyndman, R. J., & Koehler, A. B. (2006). Another look at measures of forecast accuracy. International Journal of Forecasting, 22(4), 679–688.
• Bergmeir, C., Hyndman, R. J., & Koo, B. (2018). A note on the validity of cross-validation for evaluating autoregressive time series prediction. Computational Statistics & Data Analysis, 120, 70–83.

Patterns and classification (ch. 3)

• Sakoe, H., & Chiba, S. (1978). Dynamic programming algorithm optimization for spoken word recognition. IEEE Transactions on Acoustics, Speech, and Signal Processing, 26(1), 43–49.
• Keogh, E., & Ratanamahatana, C. A. (2005). Exact indexing of dynamic time warping. Knowledge and Information Systems, 7(3), 358–386.
• Rakthanmanon, T., Campana, B., Mueen, A., Batista, G., Westover, B., Zhu, Q., Zakaria, J., & Keogh, E. (2012). Searching and mining trillions of time series subsequences under dynamic time warping. Proceedings of ACM SIGKDD, 262–270.
• Petitjean, F., Ketterlin, A., & Gançarski, P. (2011). A global averaging method for dynamic time warping, with applications to clustering. Pattern Recognition, 44(3), 678–693.
• Cuturi, M., & Blondel, M. (2017). Soft-DTW: a differentiable loss function for time-series. Proceedings of ICML, PMLR 70, 894–903.
• Cover, T., & Hart, P. (1967). Nearest neighbor pattern classification. IEEE Transactions on Information Theory, 13(1), 21–27.
• Saito, N. (1994). Local Feature Extraction and Its Applications Using a Library of Bases. PhD thesis, Yale University.
• Bagnall, A., Lines, J., Bostrom, A., Large, J., & Keogh, E. (2017). The great time series classification bake off: a review and experimental evaluation of recent algorithmic advances. Data Mining and Knowledge Discovery, 31(3), 606–660.
• Bagnall, A., Dau, H. A., Lines, J., Flynn, M., Large, J., Bostrom, A., Southam, P., & Keogh, E. (2018). The UEA multivariate time series classification archive, 2018. arXiv:1811.00075.
• Dau, H. A., Bagnall, A., Kamgar, K., Yeh, C.-C. M., Zhu, Y., Gharghabi, S., Ratanamahatana, C. A., & Keogh, E. (2019). The UCR time series archive. IEEE/CAA Journal of Automatica Sinica, 6(6), 1293–1305.
• Yeh, C.-C. M., Zhu, Y., Ulanova, L., Begum, N., Ding, Y., Dau, H. A., Silva, D. F., Mueen, A., & Keogh, E. (2016). Matrix Profile I: All pairs similarity joins for time series. Proceedings of IEEE ICDM, 1317–1322.
• Zhu, Y., Zimmerman, Z., Senobari, N. S., Yeh, C.-C. M., Funning, G., Mueen, A., Brisk, P., & Keogh, E. (2016). Matrix Profile II: Exploiting a novel algorithm and GPUs to break the one hundred million barrier for time series motifs and joins. Proceedings of IEEE ICDM, 739–748.
• Ye, L., & Keogh, E. (2009). Time series shapelets: a new primitive for data mining. Proceedings of ACM SIGKDD, 947–956.
• Grabocka, J., Schilling, N., Wistuba, M., & Schmidt-Thieme, L. (2014). Learning time-series shapelets. Proceedings of ACM SIGKDD, 392–401.

Sparse representations (ch. 4)

• Natarajan, B. K. (1995). Sparse approximate solutions to linear systems. SIAM Journal on Computing, 24(2), 227–234.
• Pati, Y. C., Rezaiifar, R., & Krishnaprasad, P. S. (1993). Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition. Proceedings of the 27th Asilomar Conference on Signals, Systems and Computers, 40–44.
• Donoho, D. L., & Elad, M. (2003). Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ1 minimization. Proceedings of the National Academy of Sciences, 100(5), 2197–2202.
• Tropp, J. A. (2004). Greed is good: algorithmic results for sparse approximation. IEEE Transactions on Information Theory, 50(10), 2231–2242.
• Beck, A., & Teboulle, M. (2009). A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences, 2(1), 183–202.
• Engan, K., Aase, S. O., & Husøy, J. H. (1999). Method of optimal directions for frame design. Proceedings of IEEE ICASSP, 2443–2446.
• Aharon, M., Elad, M., & Bruckstein, A. (2006). K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation. IEEE Transactions on Signal Processing, 54(11), 4311–4322.
• Mairal, J., Bach, F., Ponce, J., & Sapiro, G. (2010). Online learning for matrix factorization and sparse coding. Journal of Machine Learning Research, 11, 19–60.
• Dupré la Tour, T., Moreau, T., Jas, M., & Gramfort, A. (2018). Multivariate convolutional sparse coding for electromagnetic brain signals. Advances in Neural Information Processing Systems 31.

Changes and anomalies (ch. 5)

• Jackson, B., Scargle, J. D., Barnes, D., Arabhi, S., Alt, A., Gioumousis, P., Gwin, E., Sangtrakulcharoen, P., Tan, L., & Tsai, T. T. (2005). An algorithm for optimal partitioning of data on an interval. IEEE Signal Processing Letters, 12(2), 105–108.
• Killick, R., Fearnhead, P., & Eckley, I. A. (2012). Optimal detection of changepoints with a linear computational cost. Journal of the American Statistical Association, 107(500), 1590–1598.
• Fryzlewicz, P. (2014). Wild binary segmentation for multiple change-point detection. The Annals of Statistics, 42(6), 2243–2281.
• Harchaoui, Z., & Cappé, O. (2007). Retrospective multiple change-point estimation with kernels. Proceedings of the IEEE/SP Workshop on Statistical Signal Processing, 768–772.
• Arlot, S., Celisse, A., & Harchaoui, Z. (2019). A kernel multiple change-point algorithm via model selection. Journal of Machine Learning Research, 20(162), 1–56.
• Page, E. S. (1954). Continuous inspection schemes. Biometrika, 41(1/2), 100–115.
• Lorden, G. (1971). Procedures for reacting to a change in distribution. The Annals of Mathematical Statistics, 42(6), 1897–1908.
• Liu, F. T., Ting, K. M., & Zhou, Z.-H. (2008). Isolation forest. Proceedings of IEEE ICDM, 413–422.
• Kim, S., Choi, K., Choi, H.-S., Lee, B., & Yoon, S. (2022). Towards a rigorous evaluation of time-series anomaly detection. Proceedings of the AAAI Conference on Artificial Intelligence, 36.
• Wu, R., & Keogh, E. J. (2023). Current time series anomaly detection benchmarks are flawed and are creating the illusion of progress. IEEE Transactions on Knowledge and Data Engineering, 35(3), 2421–2429.

Repairing data (ch. 6)

• Rubin, D. B. (1976). Inference and missing data. Biometrika, 63(3), 581–592.
• Donoho, D. L., & Johnstone, I. M. (1994). Ideal spatial adaptation by wavelet shrinkage. Biometrika, 81(3), 425–455.
• Coifman, R. R., & Donoho, D. L. (1995). Translation-invariant de-noising. In Wavelets and Statistics, Lecture Notes in Statistics 103, Springer, 125–150.
• Rudin, L. I., Osher, S., & Fatemi, E. (1992). Nonlinear total variation based noise removal algorithms. Physica D, 60(1–4), 259–268.
• Condat, L. (2013). A direct algorithm for 1-D total variation denoising. IEEE Signal Processing Letters, 20(11), 1054–1057.
• Cai, J.-F., Candès, E. J., & Shen, Z. (2010). A singular value thresholding algorithm for matrix completion. SIAM Journal on Optimization, 20(4), 1956–1982.
• Mazumder, R., Hastie, T., & Tibshirani, R. (2010). Spectral regularization algorithms for learning large incomplete matrices. Journal of Machine Learning Research, 11, 2287–2322.
• Candès, E. J., & Recht, B. (2009). Exact matrix completion via convex optimization. Foundations of Computational Mathematics, 9(6), 717–772.

Features (ch. 7)

• Fulcher, B. D., & Jones, N. S. (2014). Highly comparative feature-based time-series classification. IEEE Transactions on Knowledge and Data Engineering, 26(12), 3026–3037.
• Lubba, C. H., Sethi, S. S., Knaute, P., Schultz, S. R., Fulcher, B. D., & Jones, N. S. (2019). catch22: CAnonical Time-series CHaracteristics. Data Mining and Knowledge Discovery, 33(6), 1821–1852.
• Christ, M., Braun, N., Neuffer, J., & Kempa-Liehr, A. W. (2018). Time series feature extraction on basis of scalable hypothesis tests (tsfresh — A Python package). Neurocomputing, 307, 72–77.
• Dempster, A., Petitjean, F., & Webb, G. I. (2020). ROCKET: exceptionally fast and accurate time series classification using random convolutional kernels. Data Mining and Knowledge Discovery, 34(5), 1454–1495.
• Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B, 58(1), 267–288.
• Kraskov, A., Stögbauer, H., & Grassberger, P. (2004). Estimating mutual information. Physical Review E, 69(6), 066138.
• Peng, H., Long, F., & Ding, C. (2005). Feature selection based on mutual information: criteria of max-dependency, max-relevance, and min-redundancy. IEEE Transactions on Pattern Analysis and Machine Intelligence, 27(8), 1226–1238.
• Meinshausen, N., & Bühlmann, P. (2010). Stability selection. Journal of the Royal Statistical Society: Series B, 72(4), 417–473.
• Ambroise, C., & McLachlan, G. J. (2002). Selection bias in gene extraction on the basis of microarray gene-expression data. Proceedings of the National Academy of Sciences, 99(10), 6562–6566.

Graphs (ch. 8)

• Granger, C. W. J. (1969). Investigating causal relations by econometric models and cross-spectral methods. Econometrica, 37(3), 424–438.
• Shuman, D. I., Narang, S. K., Frossard, P., Ortega, A., & Vandergheynst, P. (2013). The emerging field of signal processing on graphs. IEEE Signal Processing Magazine, 30(3), 83–98.
• Friedman, J., Hastie, T., & Tibshirani, R. (2008). Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3), 432–441.
• Kalofolias, V. (2016). How to learn a graph from smooth signals. Proceedings of AISTATS, PMLR 51, 920–929.
• Dong, X., Thanou, D., Frossard, P., & Vandergheynst, P. (2016). Learning Laplacian matrix in smooth graph signal representations. IEEE Transactions on Signal Processing, 64(23), 6160–6173.
• Grassi, F., Loukas, A., Perraudin, N., & Ricaud, B. (2018). A time-vertex signal processing framework: scalable processing and meaningful representations for time-series on graphs. IEEE Transactions on Signal Processing, 66(3), 817–829.

Neural networks and research directions (ch. 9–10)

• van den Oord, A., Dieleman, S., Zen, H., Simonyan, K., Vinyals, O., Graves, A., Kalchbrenner, N., Senior, A., & Kavukcuoglu, K. (2016). WaveNet: A generative model for raw audio. arXiv:1609.03499.
• Hochreiter, S., & Schmidhuber, J. (1997). Long short-term memory. Neural Computation, 9(8), 1735–1780.
• Salinas, D., Flunkert, V., Gasthaus, J., & Januschowski, T. (2020). DeepAR: Probabilistic forecasting with autoregressive recurrent networks. International Journal of Forecasting, 36(3), 1181–1191.
• Oreshkin, B. N., Carpov, D., Chapados, N., & Bengio, Y. (2020). N-BEATS: Neural basis expansion analysis for interpretable time series forecasting. International Conference on Learning Representations.
• Ismail Fawaz, H., Lucas, B., Forestier, G., Pelletier, C., Schmidt, D. F., Weber, J., Webb, G. I., Idoumghar, L., Muller, P.-A., & Petitjean, F. (2020). InceptionTime: Finding AlexNet for time series classification. Data Mining and Knowledge Discovery, 34(6), 1936–1962.
• Nie, Y., Nguyen, N. H., Sinthong, P., & Kalagnanam, J. (2023). A time series is worth 64 words: Long-term forecasting with Transformers. International Conference on Learning Representations.
• Zeng, A., Chen, M., Zhang, L., & Xu, Q. (2023). Are Transformers effective for time series forecasting? Proceedings of the AAAI Conference on Artificial Intelligence, 37.
• Ansari, A. F., Stella, L., Turkmen, C., Zhang, X., Mercado, P., Shen, H., Shchur, O., Rangapuram, S. S., Pineda Arango, S., Kapoor, S., Zschiegner, J., Maddix, D. C., Mahoney, M. W., Torkkola, K., Wilson, A. G., Bohlke-Schneider, M., & Wang, Y. (2024). Chronos: Learning the language of time series. Transactions on Machine Learning Research.
• Das, A., Kong, W., Sen, R., & Zhou, Y. (2024). A decoder-only foundation model for time-series forecasting. Proceedings of ICML, PMLR 235.
• Demšar, J. (2006). Statistical comparisons of classifiers over multiple data sets. Journal of Machine Learning Research, 7, 1–30.
• Che, Z., Purushotham, S., Cho, K., Sontag, D., & Liu, Y. (2018). Recurrent neural networks for multivariate time series with missing values. Scientific Reports, 8, 6085.
• Rubanova, Y., Chen, R. T. Q., & Duvenaud, D. (2019). Latent ordinary differential equations for irregularly-sampled time series. Advances in Neural Information Processing Systems 32.
• Tashiro, Y., Song, J., Song, Y., & Ermon, S. (2021). CSDI: Conditional score-based diffusion models for probabilistic time series imputation. Advances in Neural Information Processing Systems 34.
• Jas, M., Dupré la Tour, T., Şimşekli, U., & Gramfort, A. (2017). Learning the morphology of brain signals using alpha-stable convolutional sparse coding. Advances in Neural Information Processing Systems 30.
• Li, Y., Yu, R., Shahabi, C., & Liu, Y. (2018). Diffusion convolutional recurrent neural network: Data-driven traffic forecasting. International Conference on Learning Representations.
• Wu, Z., Pan, S., Long, G., Jiang, J., & Zhang, C. (2019). Graph WaveNet for deep spatial-temporal graph modeling. Proceedings of IJCAI, 1907–1913.