Geometric Data Analysis

Contents

• Introduction and overview
• 0. Review: linear algebra, probability and numerics
• 1. Why geometry?
• 2. Flat vector spaces
• 3. Graphs and embeddings
• 4. Geometric deep learning
• 5. Riemannian metrics and geodesics
• 6. Probability distributions
• 7. GPU programming
• 8. Complements for self-study
• 9. 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 trained a few models and wants to know which representation, distance and loss to use for vectors, graphs, surfaces and probability distributions.

Empirical first: in dimension 1000, the farthest of 1000 uniform points is only 11 % farther from a query than the nearest one (Code 1.4). Moreover, a gradient flow driven by a Gaussian kernel norm with a narrow kernel leaves every point stranded, while the same flow driven by optimal transport brings all of them to the target (Code 6.16). Both facts are geometric, and this course explains them.

The course gives one geometric viewpoint on common machine learning models. Data are vectors, points on a curved space, nodes of a graph or samples of a distribution; each choice comes with its own distances, convolutions and losses, and with its own computational cost. The course is a hub for the Santé track and for the geometry-oriented courses: Optimal Transport for Machine Learning, Geometry Processing and Geometric Deep Learning, Geometry and Shape Spaces and Graphs in Machine Learning.

What you will be able to do

• Explain the curse of dimensionality with exact computations, and say when a nearest-neighbour or kernel method can still work.
• Encode an invariance in a distance, and compute shape distances on Kendall’s sphere of triangles.
• Compare decision trees, nearest neighbours, linear models, neural networks and kernel methods by the geometry they impose on the feature space.
• Choose between a squared loss and a cross-entropy, and justify the choice by gradients and by the geometry of the simplex.
• Estimate the intrinsic dimension of a data set, measure curvature, and build graph-based embeddings (Isomap, Laplacian eigenmaps, t-SNE, UMAP).
• Write a graph convolution and a rotation-equivariant layer from scratch, and prove their equivariance.
• Compute distances, geodesics and Fréchet means on the sphere and in the Poincaré disk, and measure the error of graph shortest paths.
• Compare probability distributions with divergences, the Fisher–Rao metric, kernel norms and optimal transport, and implement the Sinkhorn algorithm.
• Relate kernel matrices, graphs, heat kernels, entropic transport plans and attention matrices to one another.
• Predict whether a computation is memory-bound or compute-bound, and write a tiled reduction that never stores an M \times N matrix.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com sums up the objective as “choosing the correct loss, distance and convolution for every problem”. It presents machine learning as the meeting point of modelling, software engineering and statistics. In that setting, geometry encodes expert knowledge: well-chosen distances between samples enforce fundamental properties while keeping a simple mathematical form.

The description stresses structured data, such as anatomical images, protein conformations and physical simulations. Its purpose is a unifying geometric perspective on common models, so that students can answer three questions. Should data be represented by vectors, graphs or histograms? How much does the choice between a Euclidean norm and a cross-entropy matter? How are kernels, graphs, manifolds, optimal transport and transformer networks related?

master-mva.com organises the course in seven lectures of three hours, five of which include a lab session. The “Thèmes abordés” section presents the class as foundational background for geometry in machine learning and lists the courses of both semesters for which it is useful. Each item is matched below.

Beyond the syllabus, this page adds a review chapter (ch. 0), a chapter of complements with longer proofs (ch. 8), a research chapter (ch. 9), graded exercises and a written-exam proposal.

Where this course fits in the programme

Official prerequisites: linear algebra and basic programming with Python (master-mva.com). The instructor listed on master-mva.com is Jean Feydy. The course carries the label of the Santé track and is taught in period 1, so only the refreshers come before it.

• Refreshers — Linear Algebra, Analysis and Optimization for the spectral theorem and gradients; Probability and Statistics for Gaussian vectors and divergences; Scientific Python, NumPy and PyTorch for the labs.
• Related courses of semester 1 — Optimal Transport for Machine Learning, Geometry Processing and Geometric Deep Learning, Topological Data Analysis, Medical Image Analysis, Deep Learning, Introduction to Statistical Learning and 3D Computer Vision.
• Courses of semester 2 that build on it — Geometry and Shape Spaces, Kernel Methods, Point Clouds and 3D Modelling, Graphs in Machine Learning, Deep Learning for Medical Imaging, Generative Models for Images, Inverse Problems and Imaging, Deformable Models and Geodesic Methods and Algorithms and Learning for Protein Science.
• Computation — Training and Deploying Large-Scale Models continues ch. 7 at the scale of large networks.
• Afterwards — the Research Internship; the Hospital Immersion of the Santé track gives a clinical context to the medical case studies.

Course outline

• Review: linear algebra, probability and numerics — inner products and the Gram trick, the spectral theorem, Gaussian vectors, softmax, cross-entropy and the stable log-sum-exp.
• Why geometry? — the curse of dimensionality, invariant distances, Kendall’s sphere of triangles, style transfer as a comparison of feature distributions.
• Flat vector spaces — trees, nearest neighbours, linear models, neural networks and kernels, read as choices of geometry on \mathbb{R}^d.
• Graphs and embeddings — Gaussian bubbles, the manifold hypothesis, intrinsic dimension, curvature, neighbourhood graphs and embeddings.
• Geometric deep learning — graph convolutions, equivariant layers, five challenges, protein docking and lung registration, learned graphs and attention.
• Riemannian metrics and geodesics — metrics, the Poincaré disk, discrete and continuous shortest paths, Fréchet means.
• Probability distributions — divergences, information geometry, kernel norms, optimal transport and the Sinkhorn algorithm; one matrix read five ways.
• GPU programming — the hardware, the roofline model, tiled reductions with an online log-sum-exp, and the end of Dennard scaling.
• Complements for self-study — the proofs left out of the sessions, Ollivier–Ricci curvature, trees in the hyperbolic plane, convergence and debiasing of Sinkhorn, the graph Laplacian as a Laplace–Beltrami operator.
• Research directions and open problems — scalable geometric computing, equivariance, protein design, unbalanced transport, transformers as interacting particles (self-study).

Session plan

21 h in class = 7 sessions of 3 h.

The volume is the official one of the Université Paris-Saclay programme table (14 h of lectures and 7 h of labs), and master-mva.com announces seven lectures of three hours each. The plan keeps one session per official lecture. Reading each session as about 2 h of lecture and 1 h of lab, which matches the 14 h + 7 h split, is a proposal, as is the content of each lab. The review chapter is not an official lecture: it is read before session 1.

Level 1 = foundations, 2 = core methods, 3 = advanced material close to research; the level never decreases from one session to the next.

Self-study, outside the 21 h (about 6.5 h in total): ch. 0 Review: linear algebra, probability and numerics, to read before session 1 (≈ 1.5 h) · ch. 8 Complements — the proofs left out of the sessions, Ollivier–Ricci curvature, hyperbolic trees, Sinkhorn convergence and debiasing, Laplacian convergence (≈ 3 h, with Ex. 17 and 18) · ch. 9 Research directions and open problems (≈ 2 h).

How to study this module

• Read ch. 0 before session 1, even if it looks familiar: the Gram trick (Proposition 0.2) and the stable log-sum-exp (Proposition 0.9) come back in ch. 3, 4, 6 and 7.
• Follow one thread through the page: a Gaussian kernel matrix is a similarity (ch. 2), a weighted graph (ch. 3), a heat kernel (ch. 5), an entropic transport plan and an attention matrix (ch. 6), and a reduction to compute fast (ch. 7).
• After each session, run the code of the chapter and change one parameter (dimension, kernel width, \varepsilon, tile size) to see what breaks.
• Do the exercises listed in the session plan before the next session, and read the proofs of ch. 8 after the session that states the result. The ★★★ exercises and the projects prepare the project report, which carries 15 of the 20 points.
• Python environment: NumPy, SciPy, scikit-learn, NetworkX and PyTorch on CPU; every code block runs in under 30 seconds. The official labs use UMAP, PyTorch Geometric, GeomStats and GeomLoss; the code here rewrites the needed parts in plain NumPy and PyTorch so that nothing else must be installed.

0. Review: linear algebra, probability and numerics

Every chapter compares data points, and every comparison is computed from inner products, eigenvalues or logarithms of sums of exponentials. This chapter states these M1 tools in the form used later, with short proofs.

0.1 Inner products and the Gram trick

All pairwise distances therefore cost one matrix product, the fastest operation on modern hardware (ch. 7).

0.2 Spectral theorem and positive matrices

A PSD similarity matrix is therefore an inner product between hidden feature vectors (ch. 2 and 3).

0.3 Gaussian vectors

0.4 Softmax, cross-entropy and the log-sum-exp

import numpy as np

rng = np.random.default_rng(0)
X, Y = rng.normal(size=(5, 3)), rng.normal(size=(4, 3))
# |x - y|^2 = |x|^2 + |y|^2 - 2 <x, y>: all distances from one matrix product (Proposition 0.2)
D2 = (X**2).sum(1)[:, None] + (Y**2).sum(1)[None, :] - 2 * X @ Y.T
print("Gram trick error:", np.abs(D2 - ((X[:, None] - Y[None]) ** 2).sum(-1)).max())

def logsumexp(s):
    """log sum_j exp(s_j), computed without overflow (Proposition 0.9)."""
    m = s.max()
    return m + np.log(np.exp(s - m).sum())

s = np.array([1000.0, 999.0, 990.0])
with np.errstate(over="ignore"):
    print("naive:", np.log(np.exp(s).sum()), "  stable:", round(logsumexp(s), 6))
p = np.exp(s - logsumexp(s))                       # softmax
q = np.array([0.5, 0.3, 0.2])
print("softmax:", np.round(p, 4), "  KL(q || p) =", round(float(np.sum(q * np.log(q / p))), 4))

1. Why geometry?

Distances between raw vectors stop being informative in high dimension, yet real data are learnable. This chapter measures the problem, then shows how a well-chosen geometry, one that encodes the invariances of the data, solves it.

1.1 The curse of dimensionality

For d = 100 and \varepsilon = 0.05, the shell holds 99.4 % of the volume.

Local methods (k-NN, kernels of width \varepsilon) average labels over such cells: with \varepsilon = 0.1 in dimension 10, there are 5^{10} \approx 10^7 cells to fill.

By Chebyshev’s inequality, most squared distances from a query then lie within a few multiples of \sqrt{1.4/d} of their mean: a few percent for d = 1000.

import numpy as np

rng = np.random.default_rng(0)
n = 1000                                   # points per data set
for d in [2, 10, 100, 1000]:
    X, q = rng.random((n, d)), rng.random(d)   # uniform cloud and a query
    dist = np.linalg.norm(X - q, axis=1)
    contrast = (dist.max() - dist.min()) / dist.min()
    print(f"d = {d:4d}:  relative contrast = {contrast:6.2f},  "
          f"mean squared distance / d = {np.mean(dist**2) / d:.3f}")

Real data escape this curse: they concentrate near low-dimensional sets (ch. 3), and their labels are invariant under known transformations, which a distance can quotient out.

1.2 Geometry as a modelling tool

Translations of images (registration), rotations of molecules (alignment before a root-mean-square deviation), relabellings of the points of a cloud (ch. 6) and reparametrizations of curves are all quotiented this way.

1.3 The sphere of triangles

A triangle is a vector z = (z_1, z_2, z_3) \in \mathbb{C}^3 of vertices, and its shape is what remains after translations, rotations and scalings z \mapsto a z + b, with a \in \mathbb{C}^* and b \in \mathbb{C}. The Helmert coordinates u = \frac{z_2 - z_1}{\sqrt 2}, \qquad v = \frac{2 z_3 - z_1 - z_2}{\sqrt 6} remove translations isometrically. Scaling is removed by imposing |u|^2 + |v|^2 = 1, and rotations act as (u, v) \mapsto e^{i\theta} (u, v).

If the three vertices are independent standard complex Gaussians, (u, v) is a standard Gaussian vector of \mathbb{C}^2. Its normalization is uniform on S^3, and the map h pushes it to the uniform law on S^2. By Archimedes’ theorem, the third coordinate is then uniform on [-1, 1].

import numpy as np

def shape_point(z):
    """Point of the unit sphere S^2 representing the triangles z (rows of 3 complex vertices)."""
    u = (z[:, 1] - z[:, 0]) / np.sqrt(2)              # Helmert coordinates: translations removed
    v = (2 * z[:, 2] - z[:, 0] - z[:, 1]) / np.sqrt(6)
    r2 = abs(u)**2 + abs(v)**2                          # scale removed by dividing by r2
    w = u * np.conj(v)                                  # rotations z -> e^{it} z leave w unchanged
    return np.c_[(abs(u)**2 - abs(v)**2), 2 * w.real, 2 * w.imag] / r2[:, None]

def kendall(z1, z2):
    """Procrustes distance rho = arccos |<z1, z2>| between centred, normalized triangles."""
    a = np.c_[(z1[1] - z1[0]) / np.sqrt(2), (2 * z1[2] - z1[0] - z1[1]) / np.sqrt(6)].ravel()
    b = np.c_[(z2[1] - z2[0]) / np.sqrt(2), (2 * z2[2] - z2[0] - z2[1]) / np.sqrt(6)].ravel()
    c = abs(np.vdot(a, b)) / (np.linalg.norm(a) * np.linalg.norm(b))
    return np.arccos(min(c, 1.0))

equi = np.array([0, 1, np.exp(1j * np.pi / 3)])
flat = np.array([0, 1, 3 + 0j])
right = np.array([0, 1, 1j])
moved = 2.5 * np.exp(0.8j) * right + (4 - 1j)            # a similar copy of `right`
P = shape_point(np.vstack([equi, np.conj(equi), flat, right, moved]))
print("sphere points (equilateral, mirror, collinear, right, moved right):")
print(np.round(P, 4))
ang = np.arccos(np.clip(P[0] @ P[3], -1, 1))
print(f"angle on S^2 = {ang:.4f},  2 * Kendall distance = {2 * kendall(equi, right):.4f}")

rng = np.random.default_rng(0)                            # Gaussian random triangles
n = 200_000
Z = rng.normal(size=(n, 3)) + 1j * rng.normal(size=(n, 3))
h3 = shape_point(Z)[:, 2]
print("quartiles of the third coordinate:", np.round(np.quantile(h3, [0.25, 0.5, 0.75]), 3))
s = np.sort(np.c_[abs(Z[:, 0] - Z[:, 1]), abs(Z[:, 1] - Z[:, 2]), abs(Z[:, 2] - Z[:, 0])]**2, axis=1)
print(f"fraction of obtuse triangles = {np.mean(s[:, 2] > s[:, 0] + s[:, 1]):.4f}")

1.4 Deep Art

Style transfer (Gatys, Ecker and Bethge, 2016) optimizes an image x so that its deep features resemble those of a content image c at one layer, and those of a style image s at several layers, through Gram matrices.

Content is compared as a vector and style as a histogram of features. The right-hand side of (b) is a kernel norm between distributions (ch. 6); Kolkin, Salavon and Shakhnarovich (2019) replaced it by optimal transport.

import torch

torch.manual_seed(0)
torch.set_num_threads(1)
n = 64
y, x = torch.meshgrid(torch.arange(n), torch.arange(n), indexing="ij")

def texture(theta):
    """Stripes of orientation theta plus a little noise: a 1 x 1 x n x n image."""
    img = torch.sin(0.6 * (x * torch.cos(theta) + y * torch.sin(theta))) + 0.3 * torch.randn(n, n)
    return img[None, None].float()

conv = torch.nn.Conv2d(1, 16, 5, padding=2, padding_mode="circular", bias=False)

def features(img):                       # 16 random filters + ReLU, flattened: 16 x n^2
    with torch.no_grad():
        return torch.relu(conv(img))[0].reshape(16, -1)

def gram(F):                             # second moments of the feature distribution
    return F @ F.T / F.shape[1]

a = texture(torch.tensor(0.4))
b = torch.roll(a, shifts=(17, 29), dims=(2, 3))                       # same texture, shifted
tiles = a.reshape(1, 1, 4, 16, 4, 16).permute(0, 1, 2, 4, 3, 5).reshape(1, 1, 16, 16, 16)
c = tiles[:, :, torch.randperm(16)].reshape(1, 1, 4, 4, 16, 16).permute(0, 1, 2, 4, 3, 5).reshape(1, 1, n, n)
d = texture(torch.tensor(1.6))                                        # other orientation

Ga = gram(features(a))
for name, img in [("shifted", b), ("tiles shuffled", c), ("other texture", d)]:
    pix = (img - a).norm() / a.norm()
    sty = (gram(features(img)) - Ga).norm() / Ga.norm()
    print(f"{name:15s}: pixel distance {pix:.3f}   Gram (style) distance {sty:.3f}")

1.5 Overview of the class

The three questions of the official objective run through this table: vectors, graphs or histograms is a choice of row; the norm against the cross-entropy is a geometry on the simplex (ch. 2, 6); kernels, graphs, manifolds, transport and transformers meet in one matrix (Section 6.5).

2. Flat vector spaces

The classical model families all work in \mathbb{R}^d, but each one imposes its own geometry. This chapter asks two questions of each model: which transformations of the features leave it unchanged, and how it decides that two points are close.

2.1 Decision trees

For a node with 5 samples of each class, G = 0.5. A split into children (4, 1) and (1, 4) gives G = 1 - 0.64 - 0.04 = 0.32 in each child, a gain of 0.18.

2.2 Nearest neighbours

At finite n, however, the distance decides which points are neighbours, and the sample size needed grows like \varepsilon^{-d} (Proposition 1.2). The rule is invariant under isometries of d, not under changes of units.

2.3 Linear models and losses

An invertible affine change of features, x \mapsto A x + c, maps the set of linear models onto itself through w \mapsto A^{-\top} w. Only the regularization \|w\|^2 sees the geometry of the features.

Geometrically, \ell_2 measures the simplex with the Euclidean metric, while \ell_{\mathrm{CE}} is a Kullback–Leibler divergence, which blows up near the faces of the simplex: confident errors cost a lot. For regression with Gaussian noise, the squared loss is the right choice, being the negative log-likelihood.

2.4 Neural networks

A ReLU network thus learns its own partition, like a tree with oblique cells. Universal approximation is the limit of this picture: Cybenko (1989) proved it for sigmoids, and Leshno, Lin, Pinkus and Schocken (1993) for every non-polynomial activation, the ReLU included. On raw vectors it treats coordinates as unrelated, and structure enters through weight sharing.

Convolutions are therefore exactly the linear layers that respect translations. Replacing translations by permutations of graph nodes or by rotations of space gives the layers of ch. 4.

2.5 Kernel methods

By Proposition 0.5, a kernel is an inner product between feature vectors \varphi(x) = k(\cdot, x). A kernel method is a linear method in the geometry of that feature space, where \|\varphi(x) - \varphi(y)\|^2 = k(x, x) + k(y, y) - 2 k(x, y) .

The width \sigma sets the geometry. When \sigma \to 0, the Gaussian kernel only sees the nearest points, as k-NN does; when \sigma is large, it behaves like a low-degree polynomial. A kernel matrix is also a weighted neighbourhood graph, the starting point of ch. 3.

import numpy as np
from sklearn.datasets import make_moons
from sklearn.model_selection import train_test_split
from sklearn.tree import DecisionTreeClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.neural_network import MLPClassifier
from sklearn.svm import SVC

X, y = make_moons(n_samples=600, noise=0.25, random_state=0)
R = np.array([[1, 1], [-1, 1]]) / np.sqrt(2)                 # rotation by 45 degrees
versions = {"original": X,
            "x2 scaled by 100": X * [1, 100],              # change of units on one feature
            "rotated by 45 deg": X @ R.T}
models = {"tree (depth 6)": DecisionTreeClassifier(max_depth=6, random_state=0),
          "15-NN": KNeighborsClassifier(15),
          "logistic": LogisticRegression(),
          "MLP 32-32": MLPClassifier((32, 32), max_iter=2000, random_state=0),
          "SVM, Gaussian": SVC(kernel="rbf", gamma=2.0)}
print(f"{'':16s}" + "".join(f"{k:>20s}" for k in versions))
for name, model in models.items():
    accs = []
    for Xv in versions.values():
        Xtr, Xte, ytr, yte = train_test_split(Xv, y, test_size=0.5, random_state=0)
        accs.append(model.fit(Xtr, ytr).score(Xte, yte))
    print(f"{name:16s}" + "".join(f"{a:20.3f}" for a in accs))

3. Graphs and embeddings

High-dimensional data are not spread uniformly: they sit near thin sets with a small intrinsic dimension and sometimes a marked curvature. This chapter measures these properties, then builds graphs that capture them and embeddings that display them.

3.1 Normal distributions are… bubbles!

The density is largest at the origin, yet almost no sample lies near it: the mass sits on a thin shell of radius \sqrt d, a bubble. Two samples are almost orthogonal, and their straight-line midpoint is atypical. In the latent space of a generative model, one interpolates along great circles of the shell rather than along chords.

import numpy as np

rng = np.random.default_rng(0)
d = 512                                            # e.g. the latent space of a generative model
X, Y = rng.normal(size=(2, 10_000, d))
nx = np.linalg.norm(X, axis=1)
print(f"|X|: mean {nx.mean():.2f}, std {nx.std():.3f}   (sqrt(d) = {np.sqrt(d):.2f})")
cos = np.sum(X * Y, axis=1) / (nx * np.linalg.norm(Y, axis=1))
print(f"cosine of two samples: mean {cos.mean():+.4f}, std {cos.std():.4f}   (1/sqrt(d) = {1/np.sqrt(d):.4f})")
mid = np.linalg.norm((X + Y) / 2, axis=1)
print(f"|midpoint|: mean {mid.mean():.2f}, i.e. {(nx.mean() - mid.mean()) / nx.std():.1f} standard deviations inside the shell")

3.2 The manifold hypothesis

The images of one object under rotations and lighting changes form such a set, with m equal to the number of degrees of freedom. On \mathcal{M}, the relevant distance follows the manifold rather than the chord, and the covering argument of Proposition 1.2 costs \varepsilon^{-m} instead of \varepsilon^{-D}. Fefferman, Mitter and Narayanan (2016) give a statistical test of the hypothesis.

3.3 Dimension

The range matters: below the noise level, data look D-dimensional. Nearest-neighbour estimators avoid choosing it.

This estimator is due to Facco, d’Errico, Rodriguez and Laio (2017). It is invariant under isometries and under scaling, and it needs no choice of radius.

3.4 Curvature

On a sphere of radius R, C(r) = 2 \pi R \sin(r / R), hence K = 1 / R^2. In the hyperbolic plane of ch. 5, C(r) = 2 \pi \sinh r, hence K = -1. Positive curvature appears with normalized embeddings and directions; negative curvature appears with hierarchies, whose number of descendants grows exponentially with depth.

For a cube, each of the 8 vertices has defect \pi/2, for a total of 4\pi = 2\pi \times 2. The angle defect is the discrete Gaussian curvature of a mesh.

3.5 Neighbourhood graphs and embeddings

Close points get close coordinates, because (f_i - f_j)^2 is paid with weight w_{ij}.

Isomap (Tenenbaum, de Silva and Langford, 2000) applies classical scaling to shortest-path distances in a k-NN graph, an estimate of geodesic distances on \mathcal{M}. When the distances are not Euclidean, B has negative eigenvalues, which measure the distortion.

The heavy tails of q and v let clusters separate in the plane. In exchange, the sizes of clusters and the gaps between them carry no metric meaning.

import numpy as np
from sklearn.datasets import load_digits, make_swiss_roll
from sklearn.decomposition import PCA
from sklearn.manifold import Isomap, SpectralEmbedding, TSNE, trustworthiness
from sklearn.model_selection import cross_val_score
from sklearn.neighbors import KNeighborsClassifier, NearestNeighbors

def two_nn(Z):
    """TwoNN estimator of the intrinsic dimension (Proposition 3.5)."""
    dist, _ = NearestNeighbors(n_neighbors=3).fit(Z).kneighbors(Z)
    mu = dist[:, 2] / dist[:, 1]                     # column 0 is the point itself
    return len(Z) / np.sum(np.log(mu))

rng = np.random.default_rng(0)
S, _ = make_swiss_roll(n_samples=4000, random_state=0)
Q, _ = np.linalg.qr(rng.normal(size=(30, 30)))       # isometric copy of the roll in R^30
print(f"TwoNN: swiss roll in R^3 {two_nn(S):.2f}, in R^30 {two_nn(np.c_[S, np.zeros((4000, 27))] @ Q.T):.2f}")

X, y = load_digits(return_X_y=True)                  # 1797 images of 8 x 8 pixels, in R^64
X = X / 16.0
print(f"TwoNN: digits in R^64 {two_nn(np.unique(X, axis=0)):.2f}")
methods = {"PCA": PCA(2),
           "Isomap (k=10)": Isomap(n_neighbors=10, n_components=2),
           "Laplacian eigenmaps": SpectralEmbedding(2, n_neighbors=10, random_state=0),
           "t-SNE": TSNE(2, perplexity=30, init="pca", random_state=0)}
for name, m in methods.items():
    Y = m.fit_transform(X)
    tw = trustworthiness(X, Y, n_neighbors=10)       # are 2D neighbours true neighbours?
    acc = cross_val_score(KNeighborsClassifier(10), Y, y, cv=5).mean()
    print(f"{name:20s} trustworthiness {tw:.3f}   10-NN accuracy in 2D {acc:.3f}")

4. Geometric deep learning

Convolutions are the linear layers that respect translations (Proposition 2.10). This chapter builds the layers that respect the symmetries of graphs and of point clouds, then tests them on two scientific problems.

4.1 Convolutions on graphs

The graph convolutional network (GCN) of Kipf and Welling (2017) keeps degree 1 with a renormalization, H' = \sigma(\hat A H W), where \hat A = \tilde D^{-1/2} (A + I) \tilde D^{-1/2} . Stacking many such layers averages features over larger and larger regions until all nodes look alike, which is called over-smoothing (Exercise 16).

4.2 Point clouds and equivariance

A point-cloud convolution sums g(x_j - x_i)\, f_j over the neighbours of x_i. For rotation invariance, the offset must be expressed in a local frame (normal and tangent directions), and tangent directions are defined only up to rotation.

4.3 Five major challenges

The official text announces five challenges without naming them. These notes propose the following five, each discussed elsewhere on the page.

1. Choosing the graph. Neighbourhoods depend on the sampling density and the scale (ch. 3).
2. Invariance. Rigid motions, local frames and orientation ambiguities (Section 4.2).
3. Discretization. A surface can be a mesh, a point cloud or a voxel grid at any resolution; predictions should not depend on that choice (ch. 6).
4. Scale. Shapes have 10^5 to 10^6 points, so neighbour searches and memory dominate (ch. 7).
5. Data. Medical data sets often hold tens of labelled shapes, so priors such as invariances must carry the information.

4.4 Case study: protein docking

Proteins interact through their surfaces, where shape and chemistry must be complementary. MaSIF (Gainza and coauthors, 2020) computes a mesh of the molecular surface, cuts geodesic patches and applies learned filters in local polar coordinates. dMaSIF (Sverrisson, Feydy, Correia and Bronstein, 2021) samples points directly on a level set of a smooth distance to the atoms. It computes chemical and curvature features there and convolves them in local frames, end to end and much faster.

A surface is the level set \mathrm{SD} = r, and points are pushed onto it by Newton-like steps along \nabla \mathrm{SD}. The construction is a log-sum-exp (Proposition 0.9), the same reduction as in ch. 6 and 7.

4.5 Case study: lung registration

Comparing inhale and exhale CT scans of a patient requires a deformation matching the two vessel trees, extracted as point clouds without correspondences; branches may be missing in either scan. Registration minimizes \mathcal{D}\big( \varphi_\sharp \mu_{\text{source}}, \mu_{\text{target}} \big) + \mathcal{R}(\varphi) over deformations \varphi. Here \mathcal{D} compares measures without correspondences (a kernel norm or a transport cost, ch. 6), and \mathcal{R} penalizes irregular deformations.

Shen, Feydy and coauthors (2021) combine robust optimal transport, which allows mass to appear or vanish where branches are missing, with learned point features and smooth deformation models.

4.6 Should we learn our graphs?

A fixed geometric graph, such as a k-NN graph in space, encodes prior knowledge at a cost O(nk). Dynamic graph CNNs (Wang and coauthors, 2019) recompute the k-NN graph in feature space at every layer. Transformers go further and learn a dense graph.

A learned dense graph is flexible but costs O(n^2) and needs data or positional encodings to recover geometry. A fixed graph is cheap and encodes geometry, but it cannot adapt to the task. Section 6.5 shows that the attention matrix is a normalized Gaussian kernel.

import networkx as nx
import torch

torch.manual_seed(0)
torch.set_num_threads(1)
G = nx.karate_club_graph()                                   # 34 members, 78 friendships
y = torch.tensor([G.nodes[i]["club"] == "Officer" for i in G]).long()
A = torch.tensor(nx.to_numpy_array(G, weight=None), dtype=torch.float32) + torch.eye(34)
dinv = A.sum(1).rsqrt()
A_hat = dinv[:, None] * A * dinv[None, :]                    # D^-1/2 (A + I) D^-1/2

class GCN(torch.nn.Module):
    def __init__(self, n_in, n_hid, n_out):
        super().__init__()
        self.W1 = torch.nn.Linear(n_in, n_hid)
        self.W2 = torch.nn.Linear(n_hid, n_out)
    def forward(self, H):                                    # two rounds of message passing
        H = torch.relu(A_hat @ self.W1(H))
        return A_hat @ self.W2(H)

X = torch.eye(34)                                            # no node features: one-hot identities
train = torch.tensor([0, 33])                                # one labelled member per club
model = GCN(34, 16, 2)
opt = torch.optim.Adam(model.parameters(), lr=0.01, weight_decay=5e-4)
for step in range(200):
    opt.zero_grad()
    loss = torch.nn.functional.cross_entropy(model(X)[train], y[train])
    loss.backward()
    opt.step()
pred = model(X).argmax(1)
print(f"accuracy on the 32 unlabelled members: {(pred == y).float()[2:].mean():.3f}")

# An E(n)-equivariant layer (Proposition 4.7): messages use only squared distances.
phi = torch.nn.Sequential(torch.nn.Linear(1, 16), torch.nn.Tanh(), torch.nn.Linear(16, 1))
def egnn_step(x):
    diff = x[:, None, :] - x[None, :, :]                     # x_i - x_j
    w = phi((diff**2).sum(-1, keepdim=True))                 # invariant weights
    return x + (diff * w).mean(1)

x = torch.randn(50, 3)
Q, _ = torch.linalg.qr(torch.randn(3, 3))                    # random orthogonal matrix
t = torch.randn(3)
with torch.no_grad():
    err = (egnn_step(x @ Q.T + t) - (egnn_step(x) @ Q.T + t)).abs().max()
print(f"equivariance error of the layer: {err:.1e}")

5. Riemannian metrics and geodesics

When data live on a curved space, straight lines leave the space and averages fall off it. This chapter replaces them by geodesics and Fréchet means, in continuous models and in their discrete approximations.

5.1 Metrics, geodesics and means

5.2 The Poincaré disk

The boundary is infinitely far: d(0, 0.9) = 2.944 and d(0, 0.99) = 5.293. Perimeters grow like e^R, like the number of nodes of a tree at depth R. Trees therefore embed in the disk with arbitrarily small distortion (Sarkar, 2011), but not in a Euclidean plane, which motivates hyperbolic embeddings of hierarchies (Nickel and Kiela, 2017).

import numpy as np

def dist(z, w):
    """Geodesic distance in the Poincare disk (Definition 5.4)."""
    return np.arccosh(1 + 2 * abs(z - w)**2 / ((1 - abs(z)**2) * (1 - abs(w)**2)))

def mobius(z, a, theta=0.0):
    """Isometry of the disk sending a to 0."""
    return np.exp(1j * theta) * (z - a) / (1 - np.conj(a) * z)

for r in [0.5, 0.9, 0.99]:
    print(f"d(0, {r}) = {dist(0, r):.4f}   log((1+r)/(1-r)) = {np.log((1 + r) / (1 - r)):.4f}")

rng = np.random.default_rng(0)
z, w = 0.5 * (rng.random(2) + 1j * rng.random(2))
a = 0.3 - 0.6j
print(f"isometry check: {dist(z, w):.6f} = {dist(mobius(z, a, 1.0), mobius(w, a, 1.0)):.6f}")

for R in [1.0, 3.0, 6.0]:                     # perimeter of the geodesic circle of radius R
    rho = np.tanh(R / 2)                      # Euclidean radius of that circle
    pts = rho * np.exp(2j * np.pi * np.arange(100_000) / 100_000)
    perim = np.sum(dist(pts, np.roll(pts, 1)))
    print(f"R = {R}: perimeter {perim:10.2f}   2 pi sinh R = {2 * np.pi * np.sinh(R):10.2f}"
          f"   2 pi R = {2 * np.pi * R:6.2f}")

5.3 Discrete and continuous models

A continuous model gives a metric tensor and geodesic equations. A discrete model gives a graph whose edges carry lengths: a k-NN graph of samples, the edges of a mesh, or a pixel grid.

The theorem needs samples spread in all directions. On a regular grid, it fails.

5.4 Shortest paths

The distance u(x) = d(x_0, x) solves the eikonal equation \|\nabla u\| = 1/F in a medium of speed F. Fast marching (Sethian, 1996) solves it in Dijkstra’s order with a consistent local update, so the error vanishes with the grid step. The heat method (Crane, Weischedel and Wardetzky, 2013) diffuses heat briefly, normalizes its gradient and integrates it, justified by Varadhan’s formula.

On a Riemannian manifold the same limit gives d_{\mathcal{M}}(x, y)^2 (Varadhan, 1967). A Gaussian kernel at small scale is a geodesic distance in disguise: kernels and manifolds meet here.

5.5 Practical applications

Hierarchies are embedded in the Poincaré disk. Directions, rotations and diffusion tensors are averaged by Fréchet means on spheres, rotation groups and positive matrices. Geodesic distances on anatomical surfaces give pose-independent descriptors. GeomStats (Miolane and coauthors, 2020) implements these spaces.

import numpy as np
from scipy.sparse import lil_matrix
from scipy.sparse.csgraph import dijkstra

def exp_map(x, v):                            # unit sphere: follow the great circle
    t = np.linalg.norm(v)
    return x if t < 1e-12 else np.cos(t) * x + np.sin(t) * v / t

def log_map(x, y):                            # inverse of exp_map at x
    c = np.clip(x @ y, -1, 1)
    u = y - c * x
    n = np.linalg.norm(u)
    return np.zeros(3) if n < 1e-12 else np.arccos(c) * u / n

rng = np.random.default_rng(0)
north = np.array([0.0, 0.0, 1.0])
pts = [exp_map(north, np.r_[0.6 * rng.normal(size=2), 0]) for _ in range(50)]
x = pts[0]
for it in range(20):                          # Riemannian gradient descent (Proposition 5.3)
    g = np.mean([log_map(x, p) for p in pts], axis=0)
    x = exp_map(x, g)
print(f"Frechet mean {np.round(x, 4)}, gradient norm {np.linalg.norm(g):.1e}")
e = np.mean(pts, axis=0)
print(f"Euclidean mean {np.round(e, 4)}, norm {np.linalg.norm(e):.4f} (off the sphere)")

n = 41                                        # grid graphs: discrete shortest paths
idx = lambda i, j: i * n + j
for name, steps in [("4-neighbour", [(0, 1), (1, 0)]),
                    ("8-neighbour", [(0, 1), (1, 0), (1, 1), (1, -1)])]:
    W = lil_matrix((n * n, n * n))
    for i in range(n):
        for j in range(n):
            for di, dj in steps:
                if 0 <= i + di < n and 0 <= j + dj < n:
                    W[idx(i, j), idx(i + di, j + dj)] = np.hypot(di, dj)
    D = dijkstra(W.tocsr(), directed=False, indices=idx(0, 0))
    for (i, j) in [(40, 0), (40, 40), (40, 20)]:
        print(f"{name}: to ({i:2d},{j:2d}) graph {D[idx(i, j)]:6.2f}, Euclidean {np.hypot(i, j):6.2f}")

6. Probability distributions

Histograms, densities, point clouds and populations are probability distributions. This chapter compares them pointwise, by blurring and by moving mass, then shows that one matrix underlies kernels, graphs, transport and attention.

6.1 Geometry in statistics

For \delta_0 and \delta_t, the values are the same for t = 0.01 and t = 100, so a loss built on them gives no gradient. Point clouds, almost surely disjoint, need the other two families.

6.2 Information geometry

Two Gaussians with the same mean and widths \sigma and 2\sigma are at distance \sqrt 2 \log 2 = 0.98, whatever \sigma.

6.3 Kernel norms

For a characteristic kernel, such as the Gaussian, the norm vanishes only if \alpha = \beta (Gretton and coauthors, 2012). The energy distance uses the conditionally positive kernel k(x, y) = -\|x - y\|. Kernel norms are cheap and smooth, but a Gaussian kernel of width \sigma hardly sees displacements much larger than \sigma.

6.4 Optimal transport

Each update makes one marginal exact; convergence is linear, with a rate that degrades as \varepsilon \to 0.

import math
import torch

torch.manual_seed(0)
torch.set_num_threads(1)
torch.set_default_dtype(torch.float64)

def sinkhorn(x, y, eps, n_iter=200):
    """Log-domain Sinkhorn for uniform weights and the cost |x - y|^2 / 2 (Algorithm 6.13)."""
    C = ((x[:, None, :] - y[None, :, :])**2).sum(-1) / 2
    la, lb = -math.log(len(x)), -math.log(len(y))           # log of the uniform weights
    f, g = torch.zeros(len(x)), torch.zeros(len(y))
    schedule = torch.logspace(math.log10(C.max()), math.log10(eps), 50).tolist() + [eps] * n_iter
    for e in schedule:                                      # epsilon-scaling, then iterations at eps
        f = -e * torch.logsumexp(lb + (g[None, :] - C) / e, dim=1)
        g = -e * torch.logsumexp(la + (f[:, None] - C) / e, dim=0)
    P = torch.exp(la + lb + (f[:, None] + g[None, :] - C) / eps)   # optimal plan
    return f.mean() + g.mean(), P, C                         # OT_eps = <a, f> + <b, g>

def sinkhorn_div(x, y, eps):
    return sinkhorn(x, y, eps)[0] - (sinkhorn(x, x, eps)[0] + sinkhorn(y, y, eps)[0]) / 2

def mmd2(x, y, kernel):
    return kernel(x, x).mean() + kernel(y, y).mean() - 2 * kernel(x, y).mean()

gauss = lambda x, y: torch.exp(-((x[:, None] - y[None])**2).sum(-1) / 2)    # sigma = 1
energy = lambda x, y: -((x[:, None] - y[None])**2).sum(-1).sqrt()          # k = -|x - y|

n = 300
x = torch.randn(n, 1)
print("   t |  KL   | Gauss MMD^2 | energy | W2 (sorted) | sqrt(2 S_eps)")
for t in [0.5, 2.0, 8.0]:
    y = torch.randn(n, 1) + t
    w2 = ((x.sort(0)[0] - y.sort(0)[0])**2).mean().sqrt()
    s = sinkhorn_div(x, y, eps=0.01)
    print(f"{t:4.1f} | {t**2 / 2:5.2f} | {mmd2(x, y, gauss):11.3f} | {mmd2(x, y, energy):6.3f} |"
          f" {w2:11.3f} | {(2 * s).clamp(min=0).sqrt():.3f}")
_, P, _ = sinkhorn(x, torch.randn(n, 1) + 2, eps=0.01)
print(f"marginal errors: {(P.sum(1) - 1 / n).abs().sum():.1e}, {(P.sum(0) - 1 / n).abs().sum():.1e}")

import math
import torch

torch.manual_seed(0)
torch.set_num_threads(1)

def sinkhorn_plan(x, y, eps, n_iter=50):
    """Optimal entropic plan for uniform weights and the cost |x - y|^2 / 2."""
    C = ((x[:, None] - y[None])**2).sum(-1) / 2
    la, lb = -math.log(len(x)), -math.log(len(y))
    f, g = torch.zeros(len(x)), torch.zeros(len(y))
    for e in torch.logspace(math.log10(C.max()), math.log10(eps), 30).tolist() + [eps] * n_iter:
        f = -e * torch.logsumexp(lb + (g[None] - C) / e, dim=1)
        g = -e * torch.logsumexp(la + (f[:, None] - C) / e, dim=0)
    return torch.exp(la + lb + (f[:, None] + g[None] - C) / eps)

def bary(P, z):                       # barycentric projection: where the plan sends each x_i
    return (P @ z) / P.sum(1, keepdim=True)

n = 200
y = torch.cat([torch.randn(n // 2, 2) * 0.15 + torch.tensor([2.0, 2.0]),     # target: two blobs
               torch.randn(n // 2, 2) * 0.15 + torch.tensor([2.0, -1.0])])
x0 = torch.rand(n, 2) * 0.8 - torch.tensor([2.0, 0.4])                      # source: a square

def flow_mmd(x, sigma=0.3, lr=2.0, steps=300):
    x = x.clone().requires_grad_(True)
    k = lambda a, b: torch.exp(-((a[:, None] - b[None])**2).sum(-1) / (2 * sigma**2))
    for _ in range(steps):
        loss = k(x, x).mean() - 2 * k(x, y).mean()          # MMD^2 up to a constant
        g, = torch.autograd.grad(loss, x)
        x = (x - lr * n * g).detach().requires_grad_(True)
    return x.detach()

def flow_sinkhorn(x, eps=0.01, lr=0.5, steps=40):
    for _ in range(steps):                                  # x_i += lr (T_xy(x_i) - T_xx(x_i))
        x = x + lr * (bary(sinkhorn_plan(x, y, eps), y) - bary(sinkhorn_plan(x, x, eps), x))
    return x

runs = [("start", x0), ("MMD, sigma=0.3", flow_mmd(x0, sigma=0.3)),
        ("MMD, sigma=2", flow_mmd(x0, sigma=2.0)), ("Sinkhorn div.", flow_sinkhorn(x0))]
for name, x in runs:
    near = torch.cdist(x, y).min(1).values                  # distance to the closest target point
    print(f"{name:14s}: median distance to target {near.median():.3f}, "
          f"points farther than 0.5: {(near > 0.5).float().mean():.0%}")

6.5 One matrix, five readings

For points x_i, y_j and a width \sigma, consider K_{ij} = \exp(-\|x_i - y_j\|^2 / 2\sigma^2).

1. It is the kernel matrix of a Gaussian kernel method (ch. 2).
2. It is the weight matrix of a neighbourhood graph (ch. 3).
3. Up to a constant factor, it is a heat kernel at time t = \sigma^2/2, whose logarithm gives geodesic distances (Proposition 5.10).
4. It is the Gibbs kernel of entropic transport with \varepsilon = \sigma^2 and cost \|x - y\|^2/2: Sinkhorn rescales its rows and columns (Proposition 6.12).
5. Normalized by rows, it is an attention matrix.

Doubly stochastic attention, which runs Sinkhorn to the end, has been studied as Sinkformers (Sander, Ablin, Blondel and Peyré, 2022).

7. GPU programming

Every method of this course reduces to matrix products, kernel sums and log-sum-exps. Their speed depends less on the number of operations than on the movement of data.

7.1 What is a GPU?

A CPU has a few cores built for latency, with large caches. A GPU has thousands of simple arithmetic units, grouped into multiprocessors that run threads in groups executing the same instruction on different data (warps of 32 threads on NVIDIA hardware).

Registers and shared memory, on each multiprocessor, are fast but hold only tens to hundreds of kilobytes. Global memory holds tens of gigabytes at a few terabytes per second on data-centre GPUs; the host link carries tens of gigabytes per second. Peak arithmetic reaches tens of TFLOP/s in single precision.

7.2 Main rules of GPU programming

This is the roofline model (Williams, Waterman and Patterson, 2009). With P_{\text{peak}} = 20 TFLOP/s and B = 2 TB/s, the ridge point is I^\star = 10 flop/byte. Adding two float32 vectors does 1 flop per 12 bytes: under 1 % of the peak.

Multiplying two n \times n matrices does 2n^3 flops for 12 n^2 bytes, an intensity n/6, past the ridge for n \geq 60. The main rules follow.

1. Keep data on the GPU: the host link is about a hundred times slower.
2. Raise the intensity: fuse operations, recompute cheap intermediates instead of storing them.
3. Expose parallelism: thousands of threads, small problems batched together.
4. Read memory contiguously; avoid branches that split a warp.
5. Use the lowest accurate precision, and stabilize reductions with the log-sum-exp.

7.3 Accelerating scientific code

Kernel sums a_i = \sum_j k(x_i, y_j)\, b_j and their log-sum-exp variants appear in ch. 2, 3, 4 and 6. Stored densely, the M \times N matrix takes 4 TB for M = N = 10^6 in float32.

For D = 3, the dense version is capped at 1.6 flop/byte, memory-bound, while blocks of 256 queries reach about 200 flop/byte. This is the design of KeOps (Charlier, Feydy, Glaunès, Collin and Durif, 2021).

Merging pairs is associative, so the reduction splits across threads. FlashAttention (Dao and coauthors, 2022) uses the same identity to compute exact attention without storing the attention matrix.

import time
import torch

torch.manual_seed(0)
torch.set_num_threads(1)
M, N, D = 4000, 6000, 3
x, y, b = torch.randn(M, D), torch.randn(N, D), torch.randn(N)

def dense_lse(x, y, b):
    """log sum_j exp(-|x_i - y_j|^2 / 2 + b_j), with the full M x N matrix in memory."""
    S = -torch.cdist(x, y)**2 / 2 + b[None, :]
    return torch.logsumexp(S, dim=1), S.numel() * 4

def tiled_lse(x, y, b, tile=512):
    """Same reduction, streamed over tiles of y with an online log-sum-exp (Proposition 7.4)."""
    m = torch.full((len(x),), -float("inf"))        # running maximum
    s = torch.zeros(len(x))                         # running sum of exp(S - m)
    for j in range(0, len(y), tile):
        S = -torch.cdist(x, y[j:j + tile])**2 / 2 + b[None, j:j + tile]
        m_new = torch.maximum(m, S.max(dim=1).values)
        s = s * torch.exp(m - m_new) + torch.exp(S - m_new[:, None]).sum(dim=1)
        m = m_new
    return m + torch.log(s), len(x) * tile * 4

for f in (dense_lse, tiled_lse):
    t0 = time.perf_counter()
    out, mem = f(x, y, b)
    print(f"{f.__name__}: largest buffer {mem / 1e6:6.1f} MB, {time.perf_counter() - t0:.3f} s")
print(f"max difference: {(dense_lse(x, y, b)[0] - tiled_lse(x, y, b)[0]).abs().max():.1e}")

# Arithmetic intensity in float32 (Proposition 7.3), with blocks of B = 256 queries:
flops, B = M * N * (3 * D + 4), 256
dense = flops / (8 * M * N)                               # the matrix is written, then read
tiled = flops / (4 * (M * D + -(-M // B) * N * (D + 1) + M))  # y read once per block of queries
print(f"arithmetic intensity: dense {dense:.2f}, tiled {tiled:.0f} flop/byte")

7.4 Is Moore’s law coming to an end?

Moore (1965) observed that the number of components per chip doubled every year, a pace later revised to about every two years. Dennard scaling (Dennard and coauthors, 1974) made smaller transistors faster at constant power density; it broke down in the mid-2000s, and clock frequencies have stayed at a few gigahertz since.

Transistors are now spent on more cores and specialized units. Leiserson and coauthors (2020) argue that future gains will come from performance engineering, algorithms and streamlined hardware. Here, speed-ups come from arithmetic intensity, not faster clocks.

8. Complements for self-study

This chapter collects the proofs left out of the sessions and five extensions of the in-class material. Each section names the chapter it complements; none is needed to follow the sessions, and all prepare the projects.

8.1 Proofs for chapters 2 and 3

Theorem 2.4 (Cover and Hart). Let X'_n be the nearest neighbour of X among n samples. On the support of the law of X, X'_n \to X almost surely, since a ball around X of any radius has positive probability and eventually contains a sample.

Given the two points, the labels Y and Y' are independent, with probabilities \eta(X) and \eta(X'_n). They disagree with probability \eta(X)\, (1 - \eta(X'_n)) + \eta(X'_n)\, (1 - \eta(X)) \longrightarrow 2\, \eta(X) (1 - \eta(X)) , by continuity of \eta; dominated convergence gives R_1. With m = \min(\eta, 1 - \eta), we have 2 \eta (1 - \eta) = 2 m (1 - m), and Jensen’s inequality for the concave map m \mapsto m (1 - m) gives R_1 \leq 2 R^\star (1 - R^\star).

Proposition 3.7 (angle defects). A face with n_f sides has angles summing to (n_f - 2)\pi, and \sum_f n_f = 2E since each edge borders two faces. All the angles of the surface therefore sum to \pi (2E - 2F), and \sum_v \delta_v = 2 \pi V - \pi (2E - 2F) = 2 \pi (V - E + F) .

Proposition 3.9 (Laplacian eigenmaps). For (a), expand (f_i - f_j)^2 and use \sum_j w_{ij} = D_{ii}. For (b), set Z = D^{1/2} Y and \tilde L = D^{-1/2} L D^{-1/2}. The problem becomes \min \operatorname{tr}(Z^\top \tilde L Z) \quad \text{subject to} \quad Z^\top Z = I, \quad Z^\top D^{1/2} \mathbf{1} = 0 . By Ky Fan’s theorem, the minimum is reached by the eigenvectors of the smallest eigenvalues of \tilde L orthogonal to D^{1/2} \mathbf{1}, which spans the kernel because the graph is connected. Going back, \tilde L z = \lambda z is L v = \lambda D v with v = D^{-1/2} z.

Proposition 3.10 (classical scaling). By Proposition 0.2, \Delta = r \mathbf{1}^\top + \mathbf{1} r^\top - 2 X X^\top with r_i = \|x_i\|^2. Since J \mathbf{1} = 0, and J X = X for centred points, -\tfrac12 J \Delta J = J X X^\top J = X X^\top . Write B = U \Lambda U^\top (Theorem 0.3): U \Lambda^{1/2} has the Gram matrix of X, hence equals X up to an orthogonal map (Proposition 0.5).

8.2 Proofs for chapters 5 and 6

Theorem 5.8 (sketch). For the upper bound, follow a minimizing geodesic and replace it by a chain of samples, each within the sampling radius of the geodesic; consecutive samples are closer than \varepsilon and the chain is barely longer than the geodesic. For the lower bound, each edge is a chord shorter than \varepsilon, much smaller than the radius of curvature, so the geodesic between its ends is at most (1 + \lambda_1') times longer. Bernstein, de Silva, Langford and Tenenbaum (2000) make the constants explicit.

Proposition 6.5 (Fisher–Rao for Gaussians). The scores are \partial_\mu \log p = \frac{x - \mu}{\sigma^2}, \qquad \partial_\sigma \log p = -\frac1\sigma + \frac{(x - \mu)^2}{\sigma^3} . Their second moments are 1/\sigma^2 and (3 - 2 + 1)/\sigma^2 = 2/\sigma^2, using \mathbb{E}(X - \mu)^4 = 3\sigma^4; the cross moment vanishes because odd moments vanish. With x = \mu / \sqrt 2, the metric becomes 2 (dx^2 + d\sigma^2)/\sigma^2, twice the hyperbolic half-plane, whose distance is \operatorname{arcosh}\Big( 1 + \frac{(x_1 - x_2)^2 + (\sigma_1 - \sigma_2)^2}{2 \sigma_1 \sigma_2} \Big) . Multiplying a metric by 2 multiplies distances by \sqrt 2.

Proposition 6.6 (natural gradient). In the coordinates \eta, the chain rule multiplies the gradient by J^\top, and the metric becomes J^\top g J. Therefore (J^\top g J)^{-1} J^\top \nabla_\theta \mathcal{L} = J^{-1} g^{-1} \nabla_\theta \mathcal{L} , which is the velocity \Delta \theta = g^{-1} \nabla_\theta \mathcal{L} expressed in the coordinates \eta.

8.3 Curvature of graphs by optimal transport

Section 3.4 measured curvature by circles and angle defects. Ollivier (2009) defines a curvature for any metric space with a random walk, using the transport distance W_1 of ch. 6.

Positive curvature means that neighbourhoods overlap, like on a sphere; negative curvature means that they drift apart, like in the hyperbolic plane.

8.4 Trees need negative curvature

For k = 8 leaves in the plane, the distortion is at least 1.83. In the hyperbolic plane, it is 1.77 for R = 2, 1.24 for R = 5 and 1.11 for R = 10. Sarkar (2011) embeds every tree in the hyperbolic plane with distortion 1 + \varepsilon.

8.5 Convergence of the Sinkhorn algorithm

For positive vectors, Hilbert’s projective metric is d_H(u, u') = \log \max_{i, j} \frac{u_i\, u'_j}{u_j\, u'_i} . It ignores positive scalings, and the coordinatewise inversion u \mapsto a / u is an isometry for it.

For \|C\|_\infty / \varepsilon = 10, the bound gives \lambda^2 = 0.99982: about 5500 iterations to divide the error by e. The bound is pessimistic, but the trend is real, and it explains the \varepsilon-scaling of Algorithm 6.13.

8.6 Why the Sinkhorn divergence is debiased

The entropic cost \mathrm{OT}_\varepsilon(\alpha, \beta) is not zero when \beta = \alpha, and its minimizer in \beta is not \alpha. When \varepsilon \to \infty, the optimal plan tends to the product \alpha \otimes \beta, so \mathrm{OT}_\varepsilon(\alpha, \beta) \to \iint c(x, y)\, d\alpha(x)\, d\beta(y) . For c = \|x - y\|^2, this is minimized by a Dirac mass at the mean of \alpha: a loss based on \mathrm{OT}_\varepsilon alone collapses the model onto the mean.

The two correction terms of S_\varepsilon remove this bias, so that S_\varepsilon(\alpha, \alpha) = 0. Exercise 15 measures the effect: for two shifted Gaussians, \sqrt{2\, \mathrm{OT}_\varepsilon} overestimates W_2 = 2.063 by 13 % at \varepsilon = 1, while \sqrt{2 S_\varepsilon} gives 2.062.

8.7 The graph Laplacian and the Laplace–Beltrami operator

Laplacian eigenmaps (Proposition 3.9) are justified by a limit theorem. Let x_1, \dots, x_n be uniform samples on a compact m-dimensional submanifold \mathcal{M}, and for a smooth f let L_t f(p) = \frac{1}{n\, t\, (4 \pi t)^{m/2}} \sum_{j = 1}^n e^{-\|p - x_j\|^2 / 4t} \big( f(p) - f(x_j) \big) . Belkin and Niyogi (2008) prove that, for a suitable t = t_n \to 0, L_t f(p) converges in probability to \frac{1}{\operatorname{vol}(\mathcal{M})}\, \Delta_{\mathcal{M}} f(p) , with \Delta_{\mathcal{M}} the Laplace–Beltrami operator taken with its non-negative sign convention.

The eigenvectors of the graph Laplacian therefore approximate the eigenfunctions of \Delta_{\mathcal{M}}, the Fourier basis of the manifold. With Varadhan’s formula (Proposition 5.10), this is the second link between Gaussian kernels and the geometry of \mathcal{M}.

9. Research directions and open problems

The tools of this course are active research topics. This chapter, for self-study, points to a few directions, with references that a project can start from.

Scalable geometric computing

Kernel sums, transport plans and attention matrices are all reductions over pairs of points. Work by Charlier, Feydy, Glaunès, Collin and Durif (2021) shows that symbolic tiled reductions run on a GPU without ever storing the M \times N matrix. FlashAttention (Dao, Fu, Ermon, Rudra and Ré, 2022) applies the same online log-sum-exp to the attention of transformers.

Open problems: reductions in high dimension, where tiles of points no longer share neighbours; and multiscale or low-rank approximations with guaranteed errors, compatible with automatic differentiation.

Geometric deep learning and equivariance

Bronstein, Bruna, Cohen and Veličković (2021) organise the field around symmetry groups and the domains on which signals live: grids, groups, graphs, geodesics and gauges. Equivariant graph networks such as those of Satorras, Hoogeboom and Welling (2021) keep the cost of message passing while respecting rigid motions.

Open problems: the trade-off between expressivity and cost when features carry higher-order tensors; and the choice of local frames on surfaces, where the tangent directions are ambiguous.

Proteins and molecules

Protein structure prediction reached experimental accuracy for many targets with AlphaFold (Jumper and coauthors, 2021). Diffusion models generate new protein backbones (Watson and coauthors, 2023), and surface-based networks such as dMaSIF (Sverrisson, Feydy, Correia and Bronstein, 2021) predict interaction sites.

Open problems: representing flexibility and conformational change, and choosing between atoms, surfaces and graphs as the primary data structure. The course Algorithms and Learning for Protein Science continues this thread.

Optimal transport beyond balanced measures

Unbalanced transport (Chizat, Peyré, Schmitzer and Vialard, 2018) lets mass be created or destroyed, which robust registration needs when structures are missing. The sample complexity of Sinkhorn divergences (Genevay, Chizat, Bach, Cuturi and Peyré, 2019) interpolates between that of kernel norms and that of transport as \varepsilon varies.

Open problems: transport between measures in high dimension, where the plug-in estimator converges slowly; and fast, differentiable solvers for the very small \varepsilon needed by sharp registrations.

Intrinsic dimension and the manifold hypothesis

Fefferman, Mitter and Narayanan (2016) give a test of the manifold hypothesis with sample complexity bounds. Pope, Zhu, Abdelkader, Goldblum and Goldstein (2021) estimate the intrinsic dimension of standard image data sets and relate it to the sample complexity of learning.

Open problems: data that mix several dimensions, noise at the scale of the neighbours (Exercise 9), and the link between intrinsic dimension and generalization of deep networks.

Hyperbolic and curved representations

Poincaré embeddings (Nickel and Kiela, 2017) and hyperbolic graph networks (Chami, Ying, Ré and Leskovec, 2019) exploit negative curvature for hierarchies. Sala, De Sa, Gu and Ré (2018) analyse the trade-off between dimension, precision and distortion.

Open problems: numerical precision near the boundary of the disk, optimization on curved spaces at scale, and spaces of mixed curvature.

Transformers as interacting particles

Geshkovski, Letrouit, Polyanskiy and Rigollet (2023) model the tokens of a transformer as particles moving under attention, a normalized Gaussian kernel (Proposition 6.17), and study their clustering. Sinkformers (Sander, Ablin, Blondel and Peyré, 2022) replace the row normalization by Sinkhorn’s.

Open problems: the role of normalization layers in these dynamics, and whether doubly stochastic attention helps on geometric data.

Exercises

The exercises are sorted by level: ★ applies the course, ★★ combines results or proves a short statement, ★★★ is open-ended. The session plan says which ones to do after each session; exercises 9 and 15 are computational.

1. ★ Shells and cells. (a) Which fraction of the volume of the unit ball lies within 0.01 of its boundary in dimensions 10, 100 and 1000? (b) How many cubes of side 0.2 are needed to cover [0, 1]^{10}, and [0, 1]^{100}?

   Solution: (a) By Proposition 1.1, 1 - 0.99^d, that is 0.096, 0.634 and 0.99996. (b) By Proposition 1.2 with \varepsilon = 0.1, at least 5^{10} \approx 9.8 \times 10^6 and 5^{100} \approx 7.9 \times 10^{69} cubes; both bounds are reached by the regular grid.

2. ★ Choosing a split. A node holds 8 positive and 4 negative samples. Split A sends (6+, 0-) left and (2+, 4-) right; split B sends (5+, 1-) left and (3+, 3-) right. Compute the Gini impurities and choose. Does the choice change if the feature is replaced by its logarithm?

   Solution: the parent has G = 1 - 4/9 - 1/9 = 0.444. Split A gives \frac12 \cdot 0 + \frac12 \cdot 0.444 = 0.222; split B gives \frac12 \cdot 0.278 + \frac12 \cdot 0.5 = 0.389. Split A wins, with a gain of 0.222 against 0.056. The logarithm is increasing, so by Proposition 2.2 the same partitions are available and the choice is unchanged.

3. ★ Two losses, two predictions. For y = 1, compute the gradients of the cross-entropy and of the squared loss with respect to the score z, at z = -3 (a mistake) and z = 3 (a correct prediction).

   Solution: by Proposition 2.7, at z = -3, p = 0.0474: the cross-entropy gives -0.9526 and the squared loss -0.0861, 11 times smaller.

   At z = 3, p = 0.9526: they give -0.0474 and -0.0043. The squared loss is weak exactly where the mistake is large.

4. ★ Is a cycle Euclidean? The 4-cycle has shortest-path distances 1 between neighbours and 2 between opposite nodes. Apply classical scaling (Proposition 3.10) and compare with the unit square, whose diagonals are \sqrt 2.

   Solution: for the cycle, B = -\frac12 J \Delta J has eigenvalues -1, 0, 2, 2; the negative eigenvalue shows that no Euclidean configuration has these distances. Indeed a node at distance 1 from two points at distance 2 is their midpoint, and both such nodes would coincide. For the square, the eigenvalues are 0, 0, 1, 1, and the two eigenvectors recover the square.

5. ★ Angle defects of the Platonic solids. Compute the angle defects of the regular tetrahedron, octahedron and icosahedron, and check Proposition 3.7.

   Solution: three, four and five equilateral triangles meet at each vertex, with defects 2\pi - \pi = \pi, 2\pi - 4\pi/3 = 2\pi/3 and 2\pi - 5\pi/3 = \pi/3. With 4, 6 and 12 vertices, the totals are all 4\pi = 2\pi \times 2, since V - E + F = 2 for these surfaces.

6. ★ Distances in the disk. Compute d(0, \frac12), d(-\frac12, \frac12) and d(\frac12, \frac{i}{2}) in the Poincaré disk, and compare the last one with the path through the origin.

   Solution: d(0, \frac12) = \log 3 = 1.099, and along the diameter distances add up: d(-\frac12, \frac12) = 2 \log 3 = 2.197. For the third pair, |z - w|^2 = \frac12 and (1 - \frac14)^2 = \frac{9}{16}, so d = \operatorname{arcosh}\big( 1 + \tfrac{16}{9} \big) = \operatorname{arcosh}(25/9) = 1.681 , shorter than 2 \log 3 through the origin: the geodesic is an arc bent towards the centre.

7. ★ Two Dirac masses. For \alpha = \delta_0 and \beta = \delta_t on the real line, compute TV, KL, the Gaussian kernel norm (width \sigma), the energy distance and W_2.

   At t = 3\sigma, compare the derivatives in t of the squared kernel norm and of the squared W_2.

   Solution: TV = 1 and KL = \infty for every t \neq 0. The squared kernel norm is 2 (1 - e^{-t^2/2\sigma^2}), the energy distance is 2|t|, and W_2 = |t|.

   At t = 3 with \sigma = 1, the derivative of the squared kernel norm is 6 e^{-4.5} = 0.067, against 2t = 6 for the squared W_2: a kernel loss barely pulls distant points.

8. ★ Rooflines. A GPU has a peak of 20 TFLOP/s and a bandwidth of 2 TB/s. Bound the throughput of a float32 vector addition and of a dense Gaussian kernel sum in dimension 3. How long does a tiled kernel sum with M = N = 10^6 take at peak?

   Solution: the ridge point is 10 flop/byte. Vector addition, at 1/12 flop/byte, is bounded by 0.17 TFLOP/s. The dense kernel sum, at 13/8 flop/byte (Proposition 7.3), is bounded by 3.25 TFLOP/s. The tiled sum is compute-bound: 13 \times 10^{12} flops take 0.65 s at peak, while the dense matrix would need 4 TB.

9. ★ Computational: noise and dimension. Sample 2000 points uniformly in the unit square, embed them in \mathbb{R}^{10}, add Gaussian noise of standard deviation s \in \{0, 0.003, 0.01, 0.03, 0.1\} in all ten directions, and estimate the dimension with the two_nn function of Code 3.12.

   Hints and expected result: with one generator seeded at 0, drawing the points first and then the noise for each level in turn, the estimates are 1.98, 4.27, 6.68, 8.28 and 9.23. The nearest-neighbour distances are about 0.01, so noise of that size already fills the ambient space at the scale TwoNN uses. Subsampling the data, which increases the neighbour distances, brings the estimate back towards 2.

10. ★★ Trees are 0-hyperbolic. A metric space is \delta-hyperbolic (Gromov) if, for all points x, y, z, w, the largest of the three sums d(x, y) + d(z, w), \quad d(x, z) + d(y, w), \quad d(x, w) + d(y, z) exceeds the second largest by at most 2\delta. Prove that a tree with its path metric is 0-hyperbolic.

    Solution: the four points span a subtree with at most two branch points. After relabelling, the paths from x and y meet at a node u, and those from z and w meet at a node v.

    Let S be the total length of the four pendant paths: S = d(x, u) + d(y, u) + d(z, v) + d(w, v) . A path from the first pair to the second crosses the segment from u to v, so d(x, y) + d(z, w) = S, \qquad d(x, z) + d(y, w) = d(x, w) + d(y, z) = S + 2\, d(u, v) . The two largest sums are equal.

11. ★★ Style losses only see second moments. Using Proposition 1.10, find two different feature distributions with the same Gram matrix. Which kernel norm of ch. 6 would separate them?

    Solution: \mu_a = \frac12 (\delta_{e_1} + \delta_{-e_1}) and \mu_b = \delta_{e_1} both have Gram matrix e_1 e_1^\top, but different means. The kernel \langle f, g \rangle^2 is not characteristic: it only sees second moments. A Gaussian kernel norm, or a transport cost, separates them, as in the optimal-transport style losses of Kolkin and coauthors (2019).

12. ★★ The simplex is a sphere. Show that the Fisher metric of the categorical distributions on n classes is \sum_i dp_i^2 / p_i, and that s_i = \sqrt{p_i} maps it to four times the round metric of the unit sphere. Deduce d_{\mathrm{FR}}(p, q) = 2 \arccos \sum_i \sqrt{p_i q_i} . Compare the distances between (0.9, 0.1) and (0.99, 0.01), and between (0.51, 0.49) and (0.6, 0.4).

    Solution: take a tangent vector dp with \sum_i dp_i = 0. The derivative of \log p_X along it is dp_X / p_X, whose second moment is the Fisher metric: \mathbb{E}\Big[ \Big( \frac{dp_X}{p_X} \Big)^2 \Big] = \sum_i p_i\, \frac{dp_i^2}{p_i^2} = \sum_i \frac{dp_i^2}{p_i} . With s_i = \sqrt{p_i}, we have dp_i = 2 s_i\, ds_i, so the metric is 4 \sum_i ds_i^2, four times the round metric. Since \sum_i s_i^2 = 1, the point s lies on the unit sphere, where the geodesic distance is \arccos \langle s, s' \rangle; the factor 4 doubles it.

    The two pairs, both moved by 0.09, are at distances 0.443 and 0.181: near the faces of the simplex, the same change costs 2.4 times more, which is why confident errors are expensive (Section 2.3).

13. ★★ Attention without positions. Prove Proposition 6.17 and Proposition 4.10. Then add sinusoidal positional encodings P_i to the inputs, h_i \mapsto h_i + P_i, and show by a two-token example that the layer is no longer permutation-equivariant.

    Hints: for the second part, take h_1 = h_2; without encodings the two outputs are equal, and swapping the tokens changes nothing. With P_1 \neq P_2, the outputs differ, and swapping the inputs does not swap the outputs, because the encodings stay attached to the positions.

14. ★★ Merging partial log-sum-exps. For pairs (m_1, s_1) and (m_2, s_2) as in Proposition 7.4, define their merge with m = \max(m_1, m_2) and s = s_1 e^{m_1 - m} + s_2 e^{m_2 - m}. Prove that the merge is associative and commutative, and deduce the depth of a parallel reduction over N blocks.

    Solution: the merge represents \log(e^{m_1} s_1 + e^{m_2} s_2), and its first component is a maximum. Both are associative and commutative, so any bracketing gives the same pair. A balanced binary tree of merges reduces N blocks in \lceil \log_2 N \rceil parallel steps.

15. ★★ Computational: the entropic bias. With the sinkhorn function of Code 6.14, take 300 samples of \mathcal{N}(0, 1) and of \mathcal{N}(2, 1), and \varepsilon \in \{1, 0.1, 0.01\}.

    Compare \sqrt{2\, \mathrm{OT}_\varepsilon(\alpha, \beta)} and \sqrt{2 S_\varepsilon(\alpha, \beta)} with the sorted W_2, and print \mathrm{OT}_\varepsilon(\alpha, \alpha).

    Hints and expected result: with seed 0, W_2 = 2.063. The biased estimate gives 2.333, 2.139 and 2.076, and \mathrm{OT}_\varepsilon(\alpha, \alpha) = 0.601, 0.163 and 0.027. The debiased one gives 2.062, 2.062 and 2.063: for a translation, the debiasing of Section 8.6 removes almost all of the entropic error.

16. ★★ Over-smoothing. Let G be connected and \hat A = \tilde D^{-1/2} (A + I) \tilde D^{-1/2}. Show that the eigenvalues of \hat A lie in (-1, 1], and that 1 is simple with eigenvector \tilde D^{1/2} \mathbf{1}.

    Deduce that \hat A^k H converges to a matrix whose rows are proportional to \sqrt{\tilde d_i} times a common vector.

    Solution: \hat A is similar to the random-walk matrix \tilde D^{-1} (A + I), which is stochastic, irreducible (connected graph) and aperiodic (self-loops). Its eigenvalue 1 is simple and the others have modulus less than 1.

    In the spectral decomposition, \hat A^k \to u u^\top with u \propto \tilde D^{1/2} \mathbf{1}, so row i of \hat A^k H tends to u_i\, u^\top H. Deep linear GCNs lose all information except the degree.

17. ★★★ Should we learn the graph? Generate point clouds of 256 points sampled on three classes of surfaces (spheres, tori, ellipsoids), at random poses and scales. Compare (a) a message-passing network on a fixed k-NN graph in space, (b) the same network with the graph recomputed in feature space at each layer, and (c) a transformer without positional encoding.

    Hints: report the accuracy against the number of training clouds, from 50 to 5000, and the training time. Expect the fixed graph to win with few data, and the learned graphs to catch up with many; test the invariance to rotations by rotating the test set. Use Proposition 4.10 to explain why (c) needs no positional encoding here.

18. ★★★ Hyperbolic versus Euclidean trees. Embed a complete binary tree of depth 6 in the Poincaré disk by placing the children of each node symmetrically at a fixed hyperbolic distance R, and in the Euclidean plane by classical scaling of the tree distances. Measure the distortion of both embeddings as functions of R.

    Hints: use Code 5.6 for distances and Proposition 8.3 for intuition. The Euclidean distortion grows with the depth, while the hyperbolic one decreases as R grows, until floating-point precision near the boundary becomes the limit (Sala and coauthors, 2018).

Projects and assessment

master-mva.com announces a project report and presentation for 15 points and a quiz for 5 points; joint projects with other classes are welcome, and the class may be validated remotely. The Université Paris-Saclay course detail describes the assessment as “Evaluation Continue non Intégrale” (partial continuous assessment). The subjects below are proposals for the project part, each sized for the 15-point report and talk. The written exam is a self-assessment proposal in the spirit of the quiz, not an official modality.

Solo projects

1. [Solo · ≈ 20 h] Paper study: Sinkhorn divergences. Reproduce the main claims of Feydy, Séjourné, Vialard, Amari, Trouvé and Peyré (2019): S_\varepsilon interpolates between transport and a kernel norm, and its gradient flows reach the target.
   • Work: extend Codes 6.14 and 6.16 to 2D shapes sampled from synthetic images; sweep \varepsilon over four decades; compare with Gaussian and energy kernel norms and with the biased \mathrm{OT}_\varepsilon.
   • Success criterion: a figure of S_\varepsilon against \varepsilon with both limits of Theorem 6.15 visible, and flows in which every point ends within one blur radius of the target.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 6 and 8.
2. [Solo · ≈ 20 h] Paper study: Poincaré embeddings. Reproduce the claim of Nickel and Kiela (2017) that hierarchies embed better in low-dimensional hyperbolic space than in Euclidean space.
   • Work: synthetic trees (balanced, random recursive, with 10^3 to 10^4 nodes); Riemannian SGD in the disk with the distance of Definition 5.4; a Euclidean baseline with the same loss.
   • Success criterion: mean rank and mean average precision of edge reconstruction in dimensions 2, 5 and 10, with the hyperbolic model ahead in dimension 2.
   • Deliverables: code, report, talk. Chapters 5 and 8.
3. [Solo · ≈ 20 h] Paper study: intrinsic dimension. Reproduce the TwoNN estimator of Facco, d’Errico, Rodriguez and Laio (2017) and its scale analysis.
   • Work: synthetic manifolds of known dimension (spheres, tori, swiss rolls, products) with noise; the scikit-learn digits; estimates on decimated subsets to scan the scale.
   • Success criterion: estimates within 10 % of the true dimension for noise-free manifolds up to dimension 10, and a curve of the estimate against the subset size that exhibits the noise scale.
   • Deliverables: code, report, talk. Chapter 3.

Group projects

1. [Group of 2–3 · ≈ 30 h each] Registration of synthetic vessel trees. Build a registration pipeline in the spirit of Shen, Feydy and coauthors (2021) on synthetic data.
   • Data: random branching curves in 3D, sampled as point clouds of 10^4 points, deformed by known smooth deformations, with 10 % of the branches removed in the target.
   • Split: one member writes the deformation model (a kernel-smoothed displacement field), one the data terms (Gaussian kernel norm, energy distance, balanced and unbalanced Sinkhorn), one the evaluation and the tiled reductions of ch. 7.
   • Success criterion: the mean error on the known displacement of the shared branches, for each data term; unbalanced transport should be the most robust to missing branches.
   • Deliverables: code, a 6–8 page report, a 15-minute talk. Chapters 4, 6 and 7.
2. [Group of 2–3 · ≈ 30 h each] Learning on synthetic molecular surfaces. Reproduce the pipeline idea of dMaSIF (Sverrisson, Feydy, Correia and Bronstein, 2021) on synthetic molecules.
   • Data: random atom clouds with two atom types; surfaces as level sets of the smooth distance of Definition 4.8; a synthetic binding-site label defined by curvature and the local atom type.
   • Split: one member samples the surfaces and computes normals and curvatures; one implements a point network and a k-NN message-passing network; one implements convolutions in local frames and tests rotation invariance.
   • Success criterion: site-prediction ROC-AUC for the three models, and the drop of each under random rotations of the test molecules (zero for the frame-based model).
   • Deliverables: code, report, talk. Chapters 3, 4 and 7.
3. [Group of 2–3 · ≈ 30 h each] A small kernel-reduction library. Implement tiled kernel sums and log-sum-exps in PyTorch, then on a free GPU tier, and measure them against the roofline.
   • Work: Gaussian and Laplace kernels, an online log-sum-exp (Proposition 7.4), automatic differentiation through the reduction; benchmarks against the dense implementation for M = N from 10^3 to 10^6.
   • Split: one member writes the reductions, one the gradients and tests, one the benchmarks and the Sinkhorn solver built on top.
   • Success criterion: agreement with the dense version to 10^{-5}, memory linear in M + N, and a roofline plot locating each kernel.
   • Deliverables: code, report, talk. Chapters 6 and 7.

Written exam (proposal)

Format: 1 hour 30 minutes, one handwritten A4 sheet allowed, no electronic devices, four exercises.

1. Question (4 points). Explain the curse of dimensionality with two exact statements, and explain why nearest-neighbour methods still work on images of handwritten digits.

   Expected: the shell fraction 1 - (1 - \varepsilon)^d and the covering bound (2\varepsilon)^{-d}, or the concentration of distances with relative spread \sqrt{1.4/d}; then the manifold hypothesis, with an intrinsic dimension near 9 for the digits (Code 3.12), and invariances encoded in the distance.

2. Question (6 points). (a) Compute d(0, 0.9) in the Poincaré disk. (b) Compute the Fisher–Rao distance between \mathcal{N}(0, 1) and \mathcal{N}(0, 4). (c) Compute W_2 between the uniform empirical measures of \{0, 1, 5\} and \{2, 3, 4\}.

   Expected: (a) \log 19 = 2.944. (b) Same mean, with \sigma_1 = 1 and \sigma_2 = 2 in Proposition 6.5: d_{\mathrm{FR}} = \sqrt 2\, \operatorname{arcosh}(1 + 1/4) = \sqrt 2 \log 2 = 0.980 .

   3. Sorting matches 0 \to 2, 1 \to 3 and 5 \to 4, so W_2^2 = (4 + 4 + 1)/3 = 3 and W_2 = 1.732.

3. Question (5 points). Prove that self-attention is permutation-equivariant, and that when all keys have the same norm, the attention weights are those of a Gaussian kernel normalized by rows.

   Expected: Proposition 4.10 (relabelling multiplies the queries, keys and values by P, and the row-wise softmax gives P \mathcal{A} P^\top), then Proposition 6.17 (expand \|q - k_j\|^2; the terms constant in j cancel in the softmax).

4. Question (5 points). A Gaussian kernel sum with M = N = 10^6 points in \mathbb{R}^3 must run on a GPU with 20 TFLOP/s and 2 TB/s. Compare a dense and a tiled implementation, and explain how to compute a stable log-sum-exp in the tiled version.

   Expected: dense storage takes 4 TB and the intensity is at most 1.6 flop/byte; tiled blocks of 256 queries reach about 200 flop/byte, compute-bound, about 0.65 s at peak (Proposition 7.3); the online pair (m, s) of Proposition 7.4, with its update and its associative merge.

Grading rubric (proposal)

References

Official references

The master-mva.com course page lists no bibliography. It links two resources of the instructor, which come first here.

• J. Feydy, course website, jeanfeydy.com/Teaching (linked from master-mva.com).
• Video recordings of the lectures, YouTube playlist (linked from master-mva.com, which notes that the audio quality improves after lecture 2).

Textbooks and monographs

• J. Feydy, Geometric Data Analysis, Beyond Convolutions, PhD thesis, Université Paris-Saclay, 2020.
• T. Hastie, R. Tibshirani and J. Friedman, The Elements of Statistical Learning, 2nd ed., Springer, 2009.
• R. Vershynin, High-Dimensional Probability: An Introduction with Applications in Data Science, Cambridge University Press, 2018.
• B. Schölkopf and A. J. Smola, Learning with Kernels, MIT Press, 2002.
• W. L. Hamilton, Graph Representation Learning, Morgan & Claypool, 2020.
• M. M. Bronstein, J. Bruna, T. Cohen and P. Veličković, “Geometric deep learning: grids, groups, graphs, geodesics, and gauges”, arXiv:2104.13478, 2021.
• J. M. Lee, Introduction to Riemannian Manifolds, 2nd ed., Graduate Texts in Mathematics 176, Springer, 2018.
• X. Pennec, S. Sommer and T. Fletcher (eds.), Riemannian Geometric Statistics in Medical Image Analysis, Academic Press, 2020.
• S. Amari, Information Geometry and Its Applications, Applied Mathematical Sciences 194, Springer, 2016.
• G. Peyré and M. Cuturi, “Computational optimal transport”, Foundations and Trends in Machine Learning 11(5–6), 2019.
• F. Santambrogio, Optimal Transport for Applied Mathematicians, Birkhäuser, 2015.
• W. W. Hwu, D. B. Kirk and I. El Hajj, Programming Massively Parallel Processors: A Hands-on Approach, 4th ed., Morgan Kaufmann, 2022.

Why geometry, flat vector spaces (ch. 1–2)

• D. G. Kendall, “Shape manifolds, Procrustean metrics, and complex projective spaces”, Bulletin of the London Mathematical Society 16(2), 1984.
• L. A. Gatys, A. S. Ecker and M. Bethge, “Image style transfer using convolutional neural networks”, CVPR, 2016.
• N. Kolkin, J. Salavon and G. Shakhnarovich, “Style transfer by relaxed optimal transport and self-similarity”, CVPR, 2019.
• T. M. Cover and P. E. Hart, “Nearest neighbor pattern classification”, IEEE Transactions on Information Theory 13(1), 1967.
• G. Cybenko, “Approximation by superpositions of a sigmoidal function”, Mathematics of Control, Signals and Systems 2(4), 1989.
• M. Leshno, V. Ya. Lin, A. Pinkus and S. Schocken, “Multilayer feedforward networks with a nonpolynomial activation function can approximate any function”, Neural Networks 6(6), 1993.
• F. Pedregosa and coauthors, “Scikit-learn: machine learning in Python”, Journal of Machine Learning Research 12, 2011.

Graphs, embeddings, dimension and curvature (ch. 3, 8)

• P. Grassberger and I. Procaccia, “Measuring the strangeness of strange attractors”, Physica D 9(1–2), 1983.
• E. Facco, M. d’Errico, A. Rodriguez and A. Laio, “Estimating the intrinsic dimension of datasets by a minimal neighborhood information”, Scientific Reports 7, 2017.
• C. Fefferman, S. Mitter and H. Narayanan, “Testing the manifold hypothesis”, Journal of the American Mathematical Society 29(4), 2016.
• J. B. Tenenbaum, V. de Silva and J. C. Langford, “A global geometric framework for nonlinear dimensionality reduction”, Science 290(5500), 2000.
• M. Belkin and P. Niyogi, “Laplacian eigenmaps for dimensionality reduction and data representation”, Neural Computation 15(6), 2003.
• M. Belkin and P. Niyogi, “Towards a theoretical foundation for Laplacian-based manifold methods”, Journal of Computer and System Sciences 74(8), 2008.
• L. van der Maaten and G. Hinton, “Visualizing data using t-SNE”, Journal of Machine Learning Research 9, 2008.
• L. McInnes, J. Healy and J. Melville, “UMAP: Uniform manifold approximation and projection for dimension reduction”, arXiv:1802.03426, 2018.
• Y. Ollivier, “Ricci curvature of Markov chains on metric spaces”, Journal of Functional Analysis 256(3), 2009.

Geometric deep learning (ch. 4)

• T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks”, ICLR, 2017.
• V. G. Satorras, E. Hoogeboom and M. Welling, “E(n) equivariant graph neural networks”, ICML, 2021.
• Y. Wang, Y. Sun, Z. Liu, S. E. Sarma, M. M. Bronstein and J. M. Solomon, “Dynamic graph CNN for learning on point clouds”, ACM Transactions on Graphics 38(5), 2019.
• A. Vaswani and coauthors, “Attention is all you need”, NeurIPS, 2017.
• P. Gainza, F. Sverrisson, F. Monti, E. Rodolà, D. Boscaini, M. M. Bronstein and B. E. Correia, “Deciphering interaction fingerprints from protein molecular surfaces using geometric deep learning”, Nature Methods 17(2), 2020.
• F. Sverrisson, J. Feydy, B. E. Correia and M. M. Bronstein, “Fast end-to-end learning on protein surfaces”, CVPR, 2021.
• Z. Shen, J. Feydy, P. Liu, A. H. Curiale, R. San José Estépar, R. San José Estépar and M. Niethammer, “Accurate point cloud registration with robust optimal transport”, NeurIPS, 2021.
• M. Fey and J. E. Lenssen, “Fast graph representation learning with PyTorch Geometric”, ICLR Workshop on Representation Learning on Graphs and Manifolds, 2019.
• W. W. Zachary, “An information flow model for conflict and fission in small groups”, Journal of Anthropological Research 33(4), 1977.

Riemannian metrics and geodesics (ch. 5, 8)

• E. W. Dijkstra, “A note on two problems in connexion with graphs”, Numerische Mathematik 1, 1959.
• M. Bernstein, V. de Silva, J. C. Langford and J. B. Tenenbaum, “Graph approximations to geodesics on embedded manifolds”, technical report, Stanford University, 2000.
• J. A. Sethian, “A fast marching level set method for monotonically advancing fronts”, Proceedings of the National Academy of Sciences 93(4), 1996.
• K. Crane, C. Weischedel and M. Wardetzky, “Geodesics in heat: a new approach to computing distance based on heat flow”, ACM Transactions on Graphics 32(5), 2013.
• S. R. S. Varadhan, “On the behavior of the fundamental solution of the heat equation with variable coefficients”, Communications on Pure and Applied Mathematics 20(2), 1967.
• M. Nickel and D. Kiela, “Poincaré embeddings for learning hierarchical representations”, NeurIPS, 2017.
• R. Sarkar, “Low distortion Delaunay embedding of trees in hyperbolic plane”, Graph Drawing (GD 2011), Lecture Notes in Computer Science 7034, Springer, 2011.
• N. Miolane and coauthors, “Geomstats: a Python package for Riemannian geometry in machine learning”, Journal of Machine Learning Research 21(223), 2020.

Probability distributions (ch. 6, 8)

• S. Amari, “Natural gradient works efficiently in learning”, Neural Computation 10(2), 1998.
• A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf and A. Smola, “A kernel two-sample test”, Journal of Machine Learning Research 13, 2012.
• M. Cuturi, “Sinkhorn distances: lightspeed computation of optimal transport”, NeurIPS, 2013.
• J. Feydy, T. Séjourné, F.-X. Vialard, S. Amari, A. Trouvé and G. Peyré, “Interpolating between optimal transport and MMD using Sinkhorn divergences”, AISTATS, 2019.
• J. Franklin and J. Lorenz, “On the scaling of multidimensional matrices”, Linear Algebra and its Applications 114–115, 1989.
• M. E. Sander, P. Ablin, M. Blondel and G. Peyré, “Sinkformers: transformers with doubly stochastic attention”, AISTATS, 2022.

GPU programming (ch. 7)

• S. Williams, A. Waterman and D. Patterson, “Roofline: an insightful visual performance model for multicore architectures”, Communications of the ACM 52(4), 2009.
• B. Charlier, J. Feydy, J. A. Glaunès, F.-D. Collin and G. Durif, “Kernel operations on the GPU, with autodiff, without memory overflows”, Journal of Machine Learning Research 22(74), 2021.
• T. Dao, D. Y. Fu, S. Ermon, A. Rudra and C. Ré, “FlashAttention: fast and memory-efficient exact attention with IO-awareness”, NeurIPS, 2022.
• A. Paszke and coauthors, “PyTorch: an imperative style, high-performance deep learning library”, NeurIPS, 2019.
• G. E. Moore, “Cramming more components onto integrated circuits”, Electronics 38(8), 1965.
• R. H. Dennard, F. H. Gaensslen, H.-N. Yu, V. L. Rideout, E. Bassous and A. R. LeBlanc, “Design of ion-implanted MOSFET’s with very small physical dimensions”, IEEE Journal of Solid-State Circuits 9(5), 1974.
• C. E. Leiserson, N. C. Thompson, J. S. Emer, B. C. Kuszmaul, B. W. Lampson, D. Sanchez and T. B. Schardl, “There’s plenty of room at the Top: what will drive computer performance after Moore’s law?”, Science 368(6495), 2020.

Research directions (ch. 9)

• J. Jumper and coauthors, “Highly accurate protein structure prediction with AlphaFold”, Nature 596, 2021.
• J. L. Watson and coauthors, “De novo design of protein structure and function with RFdiffusion”, Nature 620, 2023.
• L. Chizat, G. Peyré, B. Schmitzer and F.-X. Vialard, “Unbalanced optimal transport: dynamic and Kantorovich formulations”, Journal of Functional Analysis 274(11), 2018.
• A. Genevay, L. Chizat, F. Bach, M. Cuturi and G. Peyré, “Sample complexity of Sinkhorn divergences”, AISTATS, 2019.
• P. Pope, C. Zhu, A. Abdelkader, M. Goldblum and T. Goldstein, “The intrinsic dimension of images and its impact on learning”, ICLR, 2021.
• I. Chami, R. Ying, C. Ré and J. Leskovec, “Hyperbolic graph convolutional neural networks”, NeurIPS, 2019.
• F. Sala, C. De Sa, A. Gu and C. Ré, “Representation tradeoffs for hyperbolic embeddings”, ICML, 2018.
• B. Geshkovski, C. Letrouit, Y. Polyanskiy and P. Rigollet, “A mathematical perspective on transformers”, arXiv:2312.10794, 2023.