Refresher: Linear Algebra, Analysis and Optimization for MVA

Contents

• Introduction and overview
• 0. Review: the M1 toolkit
• 1. Eigenvalues, the spectral theorem and positive semidefinite matrices
• 2. Norms, orthogonal projections and least squares
• 3. The singular value decomposition and low-rank approximation
• 4. Matrix calculus
• 5. Hilbert spaces, orthonormal bases and Fourier analysis
• 6. Convexity, smoothness and strong convexity
• 7. Gradient descent and Newton’s method
• 8. Constraints, KKT conditions and Lagrangian duality
• 9. Stochastic gradient and proximal methods
• 10. Complements for self-study
• 11. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student entering the MVA whose linear algebra and optimisation date back a few years, or were learned without proofs, and who wants to follow the first lectures of semester 1 without stopping at every matrix identity.

Empirical first: fit a polynomial of degree 11 to 50 exact data points. Solving the normal equations returns the coefficients with a relative error of 11 %, while a QR factorisation returns them to 8 digits (ch. 2). Minimise a regularised logistic loss on the breast-cancer data of scikit-learn: after 1000 steps, gradient descent still has a gradient of norm 10^{-4}, while Newton’s method reaches 10^{-13} in 8 steps (ch. 7). Eigenvalues explain both gaps, and this page is about such explanations.

This page reviews, with proofs, the linear algebra, analysis and optimisation that the MVA course descriptions take for granted. It is independent support written for the refresher fortnight, not the official refresher course. It has two companions: Probability and Statistics and Scientific Python, NumPy and PyTorch.

The first-semester courses that lean on it most are Convex Optimization, Introduction to Statistical Learning, Theoretical Foundations of Deep Learning and Foundations of Distributed and Large-Scale Optimization.

What you will be able to do

• Diagonalise a symmetric matrix, read its spectrum and bound quadratic forms with the Rayleigh quotient.
• Recognise a positive semidefinite matrix with four equivalent tests or a Schur complement, and a positive definite one with Cholesky or Sylvester’s criterion.
• Compute operator, Frobenius and nuclear norms and condition numbers, and predict the accuracy they cost.
• Solve least-squares problems by projection, QR or SVD, and explain why the normal equations can fail.
• Compute an SVD, a pseudo-inverse and the best rank-k approximation of a matrix, and run PCA as a truncated SVD.
• Differentiate scalar functions of vectors and matrices (quadratic forms, traces, log-determinants) and check the result by automatic differentiation.
• Expand a function on an orthonormal basis of a Hilbert space, and compute Fourier series, Fourier transforms and discrete Fourier transforms with their Parseval identities.
• Prove that a function is convex, smooth or strongly convex, and compute its constants L and \mu.
• State and prove the convergence rates of gradient descent, and explain the local quadratic convergence of Newton’s method.
• Write the Lagrangian, the dual problem and the KKT conditions of a constrained problem, and solve small problems with them.
• Run stochastic gradient with a sound step-size rule, and compute proximal operators such as soft thresholding.

Official objectives and syllabus (2026–2027)

The master-mva.com planning announces refresher courses (“cours de remise à niveau”) from 14 to 24 September 2026. They are optional and without registration, and mandatory for students of the I&I track. The planning publishes neither their programme, nor their instructors, nor their volume, and the Université Paris-Saclay programme table does not list them.

This page is therefore not the official refresher. It covers the mathematical tools that the official MVA course descriptions name or assume, at the level those courses need. The 12 hours are the catalogue’s estimate of guided self-study for this page.

Since there is no official topic list, the table below plays its role. Each row is a tool. The middle column lists MVA courses whose official description (master-mva.com or the Université Paris-Saclay course detail) names a topic that relies on it, with that topic in brackets. The last column says where this page teaches it.

Every tool of the list is taught on this page. The extensions are the review of ch. 0, the self-study complements of ch. 10 (Courant–Fischer and Weyl, power iteration, the sampling theorem, acceleration, conjugate functions, a proof sketch of Slater’s theorem, and the proofs deferred from the core chapters) and the research directions of ch. 11.

Where this course fits in the programme

Official prerequisites: not displayed. The page assumes a first university course in linear algebra (matrices, determinants, eigenvalues) and in the differential calculus of several variables; ch. 0 recalls what is used.

• Companion refreshers — Probability and Statistics uses ch. 1 for covariance matrices and Gaussian vectors; Scientific Python teaches the NumPy and PyTorch of the code blocks.
• Semester 1 courses that build on it — Convex Optimization and Distributed Optimization continue ch. 6–9; Introduction to Statistical Learning and Probabilistic Graphical Models use ch. 1–4; Sub-Pixel Image Processing and Image Denoising use ch. 5; Optimal Transport uses ch. 8.
• Semester 2 courses that build on it — Kernel Methods (ch. 1, 5, 8), Graphs in Machine Learning (ch. 1, 3), Large Random Matrices (ch. 1, 3), Inverse Problems and Imaging (ch. 2, 5, 9) and Generative Models for Images (ch. 8, 9).

Course outline

• 0. Review: the M1 toolkit — linear maps, rank, inner products, trace and determinant, gradients and Taylor expansions.
• 1. Eigenvalues, the spectral theorem and positive semidefinite matrices — diagonalising symmetric matrices and testing positivity.
• 2. Norms, orthogonal projections and least squares — measuring vectors and matrices, projecting, and solving overdetermined systems stably.
• 3. The singular value decomposition and low-rank approximation — the SVD, the pseudo-inverse, the Eckart–Young theorem and PCA.
• 4. Matrix calculus — gradients of functions of matrices, the chain rule and reverse-mode differentiation.
• 5. Hilbert spaces, orthonormal bases and Fourier analysis — projection in infinite dimension, Parseval, Fourier series, the Fourier transform and the FFT.
• 6. Convexity, smoothness and strong convexity — the three properties that decide how fast optimisation algorithms converge.
• 7. Gradient descent and Newton’s method — rates for smooth and strongly convex functions, and local quadratic convergence.
• 8. Constraints, KKT conditions and Lagrangian duality — multipliers, weak and strong duality, and linear programming.
• 9. Stochastic gradient and proximal methods — optimisation with noisy gradients and with non-smooth regularisers.
• 10. Complements for self-study — min–max theorems, power iteration, sampling, acceleration, conjugates, Slater’s theorem.
• 11. Research directions and open problems — where these tools meet current research.

Session plan

12 h of guided self-study = 4 steps of 3 h. The master-mva.com planning announces the refresher fortnight (14–24 September 2026) but publishes neither its content nor its volume; the 12 h are the catalogue’s estimate for this independent page, and the split into four steps is a proposal.

Level 1 = foundations, 2 = core results and algorithms, 3 = advanced material. The level never goes down from one step to the next.

Self-study, outside the 12 h (about 6 h in total): ch. 0 Review: the M1 toolkit (≈ 1 h, before step 1, if the self-check shows gaps) · ch. 10 Complements for self-study (≈ 3 h, with Ex. 11–12) · ch. 11 Research directions and open problems (≈ 2 h, with Ex. 14).

How to study this module

• Take the diagnostic self-check first. If you answer at least eleven of the twelve questions correctly, read only the recap boxes and go to the exercises.
• Otherwise follow the four steps in order. Each chapter starts with the problem it solves and ends with a recap box. No chapter relies on a later one; ch. 10 only collects deferred proofs and extensions.
• Run every code block and change one number (a size, a step, a penalty) to see the effect predicted by the text. The code uses NumPy, SciPy, scikit-learn and PyTorch on the CPU, with synthetic data or data sets bundled with scikit-learn.
• Redo the central proofs (Theorems 1.2, 2.4, 3.4, 7.3, 7.4 and 8.4) on paper without the notes: the MVA courses take these results for granted.
• Keep ch. 10 for the week before the course that needs it; ch. 11 is for curiosity and for choosing a research project.

0. Review: the M1 toolkit

Every later chapter treats matrices as linear maps, measures vectors with inner products and differentiates functions of several variables. This chapter fixes the notation and recalls the facts that the rest of the page uses without proof.

1. Eigenvalues, the spectral theorem and positive semidefinite matrices

Covariance matrices, Hessians, kernel Gram matrices and graph Laplacians are symmetric. A symmetric matrix is diagonal in a well-chosen orthonormal basis, and this chapter characterises those with a non-negative quadratic form.

The code checks Theorem 1.2 and Proposition 1.3 on a random Gram matrix, and tests positivity with Cholesky.

import numpy as np
rng = np.random.default_rng(0)
B = rng.standard_normal((5, 3))
S = B @ B.T                          # 5 x 5 Gram matrix: PSD of rank 3
lam, U = np.linalg.eigh(S)           # ascending eigenvalues, orthonormal U
print(np.round(lam, 4))              # two eigenvalues are zero up to rounding
print(np.allclose(U @ np.diag(lam) @ U.T, S), np.allclose(U.T @ U, np.eye(5)))

# the Rayleigh quotient of random vectors stays in [lam_min, lam_max]
X = rng.standard_normal((5, 1000))
R = np.einsum("ij,ij->j", X, S @ X) / np.einsum("ij,ij->j", X, X)
print(R.min() >= lam[0] - 1e-12, R.max() <= lam[-1] + 1e-12)

# Cholesky succeeds on S + 0.1 I (positive definite) and fails on -S
L = np.linalg.cholesky(S + 0.1 * np.eye(5))
try:
    np.linalg.cholesky(-S)
except np.linalg.LinAlgError as err:
    print("-S is not positive definite:", err)

# the symmetric square root S^{1/2} = U diag(sqrt(lam)) U^T
S_half = U @ np.diag(np.sqrt(np.clip(lam, 0, None))) @ U.T
print(np.allclose(S_half @ S_half, S))
# eigenvalues [0, 0, 0.1121, 2.9218, 9.9282]; every check prints True

Three vectors of \mathbb{R}^5 have a Gram matrix of rank 3, hence two zero eigenvalues.

2. Norms, orthogonal projections and least squares

An overdetermined system Ax = b usually has no solution. Least squares replaces b by the closest point of the range of A; the condition number tells how many digits the computation loses.

Forming A^\top A squares the condition number. The QR factorisation A = QR, with orthonormal columns in Q and R upper triangular, reduces the normal equations to R x = Q^\top b and, when the residual is small, loses only about \log_{10} \kappa(A) digits.

The code fits a polynomial of degree 11 to 50 exact values, so the true coefficients are known.

import numpy as np
t = np.linspace(0, 1, 50)
A = np.vander(t, 12, increasing=True)      # columns 1, t, ..., t^11
x_true = np.ones(12)
b = A @ x_true                             # exact data: the residual is zero
print(f"cond(A) = {np.linalg.cond(A):.1e}   cond(A^T A) = {np.linalg.cond(A.T @ A):.1e}")

x_ne = np.linalg.solve(A.T @ A, A.T @ b)   # normal equations
Q, R = np.linalg.qr(A)                     # reduced QR: A = QR
x_qr = np.linalg.solve(R, Q.T @ b)
x_sv = np.linalg.lstsq(A, b, rcond=None)[0]  # SVD-based solver
for name, x in [("normal equations", x_ne), ("QR", x_qr), ("lstsq (SVD)", x_sv)]:
    err = np.linalg.norm(x - x_true) / np.linalg.norm(x_true)
    print(f"{name:17s} relative error {err:.1e}")
# cond(A) = 1.2e+08, cond(A^T A) = 1.4e+16
# relative errors: normal equations 1.1e-01, QR 7.5e-09, lstsq (SVD) 3.1e-09

With \kappa(A) \approx 10^8, QR and the SVD keep eight digits; the normal equations, with \kappa(A)^2 \approx 10^{16}, keep none.

3. The singular value decomposition and low-rank approximation

The singular value decomposition (SVD) diagonalises any rectangular matrix with two orthonormal bases. Its truncation is the best low-rank approximation: the mathematics of PCA, compression and low-rank fine-tuning.

The range of A is spanned by u_1, \dots, u_r, and its kernel is orthogonal to v_1, \dots, v_r. Moreover, \|A\|_2 = \sigma_1, \qquad \|A\|_F^2 = \sum_i \sigma_i^2, \qquad \kappa(A) = \frac{\sigma_1}{\sigma_n} \ \text{ when } r = n .

The code runs PCA on the 1797 digits of scikit-learn (64 pixels each) and checks Theorem 3.4.

import numpy as np
from sklearn.datasets import load_digits

X = load_digits().data                     # 1797 images of 8 x 8 = 64 pixels
Xc = X - X.mean(axis=0)                    # centre every pixel
U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
for k in (2, 10, 30):
    Xk = (U[:, :k] * s[:k]) @ Vt[:k]       # truncated SVD = best rank-k approximation
    ok_F = np.isclose(np.linalg.norm(Xc - Xk), np.sqrt((s[k:] ** 2).sum()))
    ok_2 = np.isclose(np.linalg.norm(Xc - Xk, 2), s[k])
    print(k, ok_F, ok_2, f"explained variance {(s[:k] ** 2).sum() / (s ** 2).sum():.3f}")

# a random 10-dimensional subspace does much worse than the top 10 singular vectors
rng = np.random.default_rng(0)
P = np.linalg.qr(rng.standard_normal((64, 10)))[0]
print(f"random subspace error {np.linalg.norm(Xc - Xc @ P @ P.T):.0f}",
      f"vs Eckart-Young {np.sqrt((s[10:] ** 2).sum()):.0f}")
# explained variance 0.285 (k = 2), 0.738 (k = 10), 0.959 (k = 30); errors 1329 vs 752

Ten components out of 64 keep 74 % of the variance; a random subspace of the same dimension does far worse.

4. Matrix calculus

Training a model means differentiating a scalar loss with respect to vectors and matrices of parameters. One method, the differential, turns every such computation into a few lines of algebra.

Automatic differentiation checks such formulas in a few lines.

import torch
torch.set_num_threads(1)
torch.manual_seed(0)
f64 = dict(dtype=torch.float64)

n = 4
M = torch.randn(n, n, **f64)
X = (M @ M.T + n * torch.eye(n, **f64)).requires_grad_()   # X is positive definite
torch.logdet(X).backward()
print(torch.allclose(X.grad, torch.linalg.inv(X).T))       # gradient of log det X is X^{-T}

A, B, C = torch.randn(3, n, **f64), torch.randn(n, 2, **f64), torch.randn(3, 2, **f64)
Y = torch.randn(n, n, **f64, requires_grad=True)
((A @ Y @ B - C) ** 2).sum().backward()                    # squared Frobenius norm
print(torch.allclose(Y.grad, 2 * A.T @ (A @ Y.detach() @ B - C) @ B.T))

Both lines print True.

5. Hilbert spaces, orthonormal bases and Fourier analysis

Signals, images and the functions learned by kernel methods live in spaces of infinite dimension. Hilbert spaces keep the geometry of \mathbb{R}^n, and the Fourier basis is their most useful orthonormal basis.

The code checks Proposition 5.10 and Example 5.7.

import numpy as np
rng = np.random.default_rng(0)
N = 256
x, y = rng.standard_normal(N), rng.standard_normal(N)
X = np.fft.fft(x)
print(np.isclose((np.abs(X) ** 2).sum() / N, (x ** 2).sum()))        # Parseval for the DFT
circ = np.array([sum(x[m] * y[(k - m) % N] for m in range(N)) for k in range(N)])
print(np.allclose(np.fft.ifft(X * np.fft.fft(y)).real, circ))       # convolution theorem

# partial Fourier sums of the sawtooth f(t) = t on (-pi, pi): b_n = 2 (-1)^(n+1) / n
t = np.linspace(-np.pi, np.pi, 20001)
for K in (10, 100, 1000):
    S = sum(2 * (-1) ** (n + 1) / n * np.sin(n * t) for n in range(1, K + 1))
    rms = np.sqrt(np.mean((S - t) ** 2))
    print(f"K = {K:4d}   max S_K = {S.max():.3f}   L2 error = {rms:.3f}")
# K = 10: max 3.414, L2 error 0.437;  K = 100: 3.673, 0.143;  K = 1000: 3.701, 0.050

The L^2 error tends to 0, but the maximum tends to 3.70, not \pi: the excess is 9 % of the jump 2\pi.

6. Convexity, smoothness and strong convexity

Three properties of a function decide how fast the algorithms of the next chapters converge: convexity, smoothness and strong convexity.

Non-negative sums, affine substitutions and suprema preserve convexity: norms, the largest eigenvalue of a symmetric matrix and the losses of ch. 4 are convex.

7. Gradient descent and Newton’s method

Gradient descent uses first derivatives, and its speed is set by the condition number; Newton’s method also uses the Hessian and converges quadratically near the solution.

No convexity is needed, but only small gradients are guaranteed, not a minimiser: this is the setting of deep learning.

Nesterov’s accelerated gradient improves the two rates to 1/K^2 and (1 - 1/\sqrt\kappa)^k, which is optimal for first-order methods (ch. 10).

Far from x^*, Newton’s method needs a line search. The code compares both methods on the logistic regression of Example 6.6.

import numpy as np
from sklearn.datasets import load_breast_cancer

X, y = load_breast_cancer(return_X_y=True)          # 569 tumours, 30 features
X = (X - X.mean(0)) / X.std(0)
X = np.c_[np.ones(len(X)), X]                        # intercept column
n, d = X.shape
lam = 1e-2                                           # ridge term: mu = lam

def grad_hess(w):
    p = 1 / (1 + np.exp(-X @ w))
    g = X.T @ (p - y) / n + lam * w
    H = (X.T * (p * (1 - p))) @ X / n + lam * np.eye(d)
    return g, H

L = np.linalg.norm(X, 2) ** 2 / (4 * n) + lam        # smoothness constant
print(f"L = {L:.2f}, mu = {lam}, kappa <= {L / lam:.0f}")
w = np.zeros(d)
for k in range(1, 1001):
    w -= grad_hess(w)[0] / L                         # gradient descent, step 1/L
    if k in (10, 100, 1000):
        print(f"GD     k = {k:4d}   |grad| = {np.linalg.norm(grad_hess(w)[0]):.1e}")
w = np.zeros(d)
for k in range(1, 9):
    g, H = grad_hess(w)
    w -= np.linalg.solve(H, g)                       # pure Newton step
    print(f"Newton k = {k:4d}   |grad| = {np.linalg.norm(grad_hess(w)[0]):.1e}")
# GD: 1.3e-01, 1.4e-02, 1.2e-04 after 10, 100, 1000 steps
# Newton: 4.0e-01, 1.6e-01, 5.9e-02, 1.8e-02, 2.8e-03, 9.5e-05, 1.1e-07, 1.4e-13

Gradient descent gains about one digit per 450 steps. Newton’s method takes four slow steps, then enters the quadratic phase of Theorem 7.7.

8. Constraints, KKT conditions and Lagrangian duality

Multipliers put a price on each constraint, and the best prices solve a second problem, the dual. Consider \text{(P)} \qquad \min_x\ f_0(x) \quad \text{subject to} \quad f_i(x) \le 0 \ \ (i \le m), \qquad h_j(x) = 0 \ \ (j \le p), with optimal value p^*.

The code solves a small transport problem and reads the dual variables returned by the solver.

import numpy as np
from scipy.optimize import linprog

a, b = np.array([0.5, 0.3, 0.2]), np.array([0.4, 0.4, 0.2])   # two histograms
xs, ys = np.array([0.0, 1.0, 2.0]), np.array([0.5, 1.5, 3.0])
C = (xs[:, None] - ys[None, :]) ** 2                            # transport costs
# primal: minimise <C, P> subject to P 1 = a, P^T 1 = b, P >= 0
A_eq = np.vstack([np.kron(np.eye(3), np.ones(3)), np.kron(np.ones(3), np.eye(3))])
res = linprog(C.ravel(), A_eq=A_eq, b_eq=np.r_[a, b], bounds=(0, None), method="highs")
f, g = res.eqlin.marginals[:3], res.eqlin.marginals[3:]          # dual variables
P = res.x.reshape(3, 3)
print(f"primal value {res.fun:.4f}   dual value {a @ f + b @ g:.4f}")
print("dual feasible:", np.all(f[:, None] + g[None, :] <= C + 1e-9))
print("complementary slackness:", np.allclose(P * (C - f[:, None] - g[None, :]), 0))
print(np.round(P, 2) + 0.0)
# primal value 0.6000, dual value 0.6000; both checks print True

Mass moves only along pairs where the dual constraint is tight, f_i + g_j = C_{ij}.

9. Stochastic gradient and proximal methods

Learning objectives average many examples and often add a non-smooth regulariser such as the \ell_1 norm. Stochastic gradient handles the first, proximal operators the second.

On the whole space a bounded second moment is impossible; a projection onto a ball containing x^* enforces it without changing the proof.

import numpy as np

rng = np.random.default_rng(1)
n, d = 2000, 10
A = rng.standard_normal((n, d))
b = A @ rng.standard_normal(d) + 0.5 * rng.standard_normal(n)
x_hat = np.linalg.lstsq(A, b, rcond=None)[0]      # minimiser of f = (1/2n) |Ax - b|^2
mu = np.linalg.eigvalsh(A.T @ A / n).min()        # strong convexity constant of f

def sgd(step, iters=50_000, seed=2):
    r, x = np.random.default_rng(seed), np.zeros(d)
    for k in range(iters):
        i = r.integers(n)                         # one example, drawn uniformly
        g = (A[i] @ x - b[i]) * A[i]              # unbiased estimate of grad f(x)
        x -= step(k) * g
    return np.sum((x - x_hat) ** 2)

print(f"mu = {mu:.2f}")
for name, step in [("constant 0.02 ", lambda k: 0.02),
                   ("constant 0.002", lambda k: 0.002),
                   ("1/(mu(k+k0))  ", lambda k: 1 / (mu * (k + 100)))]:
    print(f"{name}  |x_K - x_hat|^2 = {sgd(step):.1e}")
# mu = 0.90;  constant 0.02: 1.3e-02;  constant 0.002: 2.8e-03;  decreasing: 3.3e-05

With G^2 \approx 0.25 \times 10, the noise variance times \mathbb{E}\|a_i\|^2, the floor \eta G^2 / (2\mu) is 2.8 \times 10^{-2} \ \text{for } \eta = 0.02, \qquad 2.8 \times 10^{-3} \ \text{for } \eta = 0.002 . Single runs scatter around it, within a factor of 3. The shift k_0 stabilises the first steps.

The code recovers a sparse signal from fewer measurements than unknowns with ISTA.

import numpy as np
from sklearn.linear_model import Lasso

rng = np.random.default_rng(0)
n, d = 100, 200                                  # more unknowns than equations
A = rng.standard_normal((n, d)) / np.sqrt(n)
x_true = np.zeros(d)
x_true[:5] = [3, -2, 2, -1.5, 1]                 # a 5-sparse signal
b = A @ x_true + 0.05 * rng.standard_normal(n)
lam = 0.25

def soft(v, t):                                  # prox of t * |.|_1
    return np.sign(v) * np.maximum(np.abs(v) - t, 0)

def F(x):
    return 0.5 * np.sum((A @ x - b) ** 2) + lam * np.abs(x).sum()

L = np.linalg.norm(A, 2) ** 2                    # Lipschitz constant of the smooth part
x = np.zeros(d)
for k in range(1, 3001):
    x = soft(x - A.T @ (A @ x - b) / L, lam / L) # ISTA step
    if k in (10, 100, 3000):
        print(f"k = {k:4d}   F = {F(x):.6f}")
ref = Lasso(alpha=lam / n, fit_intercept=False, tol=1e-12, max_iter=100000).fit(A, b).coef_
print(f"scikit-learn F = {F(ref):.6f}")
print("support:", np.flatnonzero(np.abs(x) > 1e-6), "values:", np.round(x[:5], 2))
# F = 2.965473, 2.343866, 2.343866 after 10, 100, 3000 steps; scikit-learn 2.343866
# support [0 1 2 3 4], values [2.83 -1.77 1.71 -1.26 0.69]

ISTA recovers the support exactly, with values shrunk towards 0: the bias of the \ell_1 penalty.

10. Complements for self-study

This chapter collects the proofs deferred from the core chapters and further tools: min–max theorems for eigenvalues, the Schur complement, power iteration, the Riesz and sampling theorems, acceleration and convex conjugates. Each section is independent; read it before the MVA course that needs it.

Min–max theorems and the Frobenius case of Eckart–Young

Proof of the Frobenius case of Theorem 3.4. Let B have rank at most k.

Weyl’s inequality with X = A - B, Y = B and j = k + 1 gives \sigma_{i+k}(A) \le \sigma_i(A - B), because \sigma_{k+1}(B) = 0.

Summing the squares over i, \|A - B\|_F^2 = \sum_{i \ge 1} \sigma_i(A - B)^2 \ \ge\ \sum_{i \ge 1} \sigma_{i+k}(A)^2 = \|A - A_k\|_F^2 . ∎

The Schur complement

The Schur complement is the conditional covariance of a Gaussian vector, and it turns constraints such as x^\top A^{-1} x \le t into the linear matrix inequalities of semidefinite programming.

Power iteration

Power iteration computes PageRank and the top singular vectors of huge matrices. Its block version, with a random starting block, is the core of the randomized SVD of Halko, Martinsson and Tropp (2011); Krylov methods such as Lanczos converge faster by keeping all the iterates.

Deferred proofs of chapter 5, and the Riesz theorem

Full proof of Theorem 5.2. Let d be the distance from x to C, and y_k \in C with \|x - y_k\| \to d.

1. The midpoint of y_k and y_l lies in C, so it is at distance at least d from x. The parallelogram law, applied to the two vectors from x, gives \|y_k - y_l\|^2 = 2\|x - y_k\|^2 + 2\|x - y_l\|^2 - 4\Big\|x - \tfrac{y_k + y_l}{2}\Big\|^2 \le 2\|x - y_k\|^2 + 2\|x - y_l\|^2 - 4d^2 .
2. The right side tends to 0, so (y_k) is a Cauchy sequence. Its limit p exists by completeness, lies in C because C is closed, and satisfies \|x - p\| = d.
3. If p and p' are two minimisers, the alternating sequence p, p', p, p', \dots is minimising, hence Cauchy, so p = p'.
4. For y \in C and t \in (0, 1], the point p + t(y - p) lies in C. Expanding \|x - p - t(y - p)\|^2 \ge \|x - p\|^2 gives -2t \operatorname{Re} \langle x - p,\ y - p \rangle + t^2 \|y - p\|^2 \ge 0 ; divide by t and let t \to 0. ∎

In a Hilbert space of functions where every evaluation f \mapsto f(x) is continuous, Riesz gives a K_x with f(x) = \langle f, K_x \rangle; the function K(x, y) = \langle K_y, K_x \rangle is the reproducing kernel of Kernel Methods.

The sampling theorem

Sampling below the rate 2B folds high frequencies onto low ones: this is aliasing, the trade-off studied in Sub-Pixel Image Processing. Mallat (2009) treats sampling, filtering and wavelet bases in depth.

Acceleration and lower bounds

These rates cannot be improved by any method that only combines gradients. Nesterov (2018, section 2.1.2) builds, for every such method and every k below half the dimension, an L-smooth convex quadratic with f(x_k) - f^* \ \ge\ \frac{3 L\, \|x_0 - x^*\|^2}{32\, (k + 1)^2} . Exercise 12 measures the gain: on a quadratic with \kappa = 1000, acceleration needs about 19 times fewer iterations than gradient descent.

Conjugate functions and proximal operators

For g = \lambda \|\cdot\|_1, the conjugate is the indicator of the box \{\|y\|_\infty \le \lambda\}, whose prox is the clipping of each coordinate to [-\lambda, \lambda]. Moreau’s identity then says that soft thresholding is v minus its clipping, which is the formula of Proposition 9.5.

Slater’s theorem

Proof sketch of Theorem 8.3, for inequality constraints only. Consider the set of achievable constraint and objective levels, \mathcal{A} = \{ (u, t) \in \mathbb{R}^m \times \mathbb{R} : f_i(x) \le u_i \text{ for all } i \text{ and } f_0(x) \le t \text{ for some } x \}.

1. \mathcal{A} is convex, because the f_i are convex, and it does not meet the half-line of points (0, s) with s < p^*.
2. A separating hyperplane gives (\lambda, \lambda_0) \neq 0 with \lambda^\top u + \lambda_0 t \ge \lambda_0 p^* on \mathcal{A}. Since the set is closed under increasing u and t, both multipliers are non-negative.
3. If \lambda_0 = 0, then \lambda^\top u \ge 0 on \mathcal{A}, which fails at the strictly feasible point, where every u_i < 0. So \lambda_0 > 0.
4. Dividing by \lambda_0: for every x, f_0(x) + \sum_i (\lambda_i / \lambda_0) f_i(x) \ge p^*. So the dual function at \lambda / \lambda_0 is at least p^*, and weak duality gives equality. ∎

Proof of Theorem 7.7

Write e_k = x_k - x^* and H_k = \nabla^2 f(x_k). The Lipschitz bound gives H_k \succeq (\mu - M \|e_k\|) I \succeq \frac{\mu}{2} I while \|e_k\| \le \mu / (2M). Since \nabla f(x^*) = 0, Taylor’s formula with integral remainder gives e_{k+1} = e_k - H_k^{-1} \nabla f(x_k) = H_k^{-1} \int_0^1 \big( H_k - \nabla^2 f(x^* + t e_k) \big)\, e_k\, dt . The integrand has norm at most M (1 - t) \|e_k\|^2, hence \|e_{k+1}\| \le \frac{2}{\mu} \cdot \frac{M}{2} \|e_k\|^2 = \frac{M}{\mu} \|e_k\|^2 \le \frac12 \|e_k\| , so the iterates stay in the region and the bound holds at every step. ∎

11. Research directions and open problems

The tools of this page are old, but their use in learning raises open questions. This chapter, for self-study, points to a few active lines of work, with entry papers.

Randomized numerical linear algebra

Halko, Martinsson and Tropp (2011) showed that a random projection followed by a small exact SVD approximates the top singular subspace of a large matrix, with error close to \sigma_{k+1} after a few power iterations. Martinsson and Tropp (2020) survey the field: sketching, randomized least squares, Nyström approximations of kernel matrices. Open questions include the right sketch size for a given accuracy, and the behaviour of these algorithms in the low-precision arithmetic of GPUs.

Low rank in learning

Candès and Recht (2009) proved that an incoherent low-rank matrix can be recovered exactly, with high probability, from few random entries by minimising the nuclear norm (Definition 3.2), a convex problem. Practitioners rather factorise X = UV^\top and run gradient descent, a non-convex problem. Burer and Monteiro (2003) introduced this factorisation for semidefinite programs. Ge, Lee and Ma (2016) showed that, for positive semidefinite matrix completion under incoherence assumptions, the factorised objective has no spurious local minimum.

Gunasekar and coauthors (2017) conjectured that gradient descent on a factorised model implicitly minimises the nuclear norm; Arora, Cohen, Hu and Luo (2019) gave evidence that deep factorisations rather favour low rank in a way no norm captures. Low-rank adaptation (Hu and coauthors, 2022) uses the same structure to fine-tune large language models, as taught in Foundations of LLMs and LLMs for Code and Proof.

The implicit bias of gradient descent

On linearly separable data, the logistic loss has no minimiser, yet Soudry and coauthors (2018) proved that gradient descent converges in direction to the maximum-margin separator, at a logarithmic rate. Which solution an optimiser selects, among the many that fit the data, is now a central question of the theory of deep learning taught in Theoretical Foundations of Deep Learning.

Beyond the classical step-size theory

The descent lemma guarantees progress only for steps below 2/L. Cohen and coauthors (2021) observed that full-batch gradient descent on neural networks drives the largest Hessian eigenvalue up to about 2/\eta and keeps training at this “edge of stability”, where the loss decreases non-monotonically. Explaining this regime is open. A related line weakens convexity: with smoothness, the Polyak–Łojasiewicz inequality of Proposition 6.5 gives linear convergence without convexity (Polyak, 1963; Karimi, Nutini and Schmidt, 2016), and holds for some over-parameterised models.

Worst-case analysis by computer

Drori and Teboulle (2014) showed that the worst case of a first-order method over all smooth convex functions can be bounded by solving a semidefinite program. Taylor, Hendrickx and Glineur (2017) made this exact with interpolation conditions for smooth strongly convex functions. The approach has produced tight constants and new methods, and it turns the design of optimisation algorithms into a computational question.

Stochastic and adaptive methods

For finite sums, variance reduction recovers linear rates with stochastic steps: SAG (Schmidt, Le Roux and Bach, 2017), SVRG (Johnson and Zhang, 2013) and SAGA (Defazio, Bach and Lacoste-Julien, 2014). For least squares, Bach and Moulines (2013) showed that averaged SGD with a constant step reaches the optimal O(1/n) rate without strong convexity.

Adaptive methods such as Adam (Kingma and Ba, 2015) and matrix preconditioners such as Shampoo (Gupta, Koren and Singer, 2018) dominate the training of large models, but their theory lags behind practice. The optimisers used for language models, discussed in Foundations of LLMs and Training and Deploying Large-Scale Models, are a fast-moving research topic.

Proximal algorithms with learned denoisers

Plug-and-play methods (Venkatakrishnan, Bouman and Wohlberg, 2013) replace the prox of a regulariser in a splitting algorithm by an image denoiser, possibly a neural network. Convergence then needs the denoiser to behave like a prox: Ryu and coauthors (2019) imposed Lipschitz constraints, and Hurault, Leclaire and Papadakis (2022) built denoisers that are exact gradient steps. These methods are studied in Generative Models for Images and Deep Learning for Image Restoration.

Exercises

The exercises are sorted by level: ★ applies the course, ★★ combines results or proves a short one, ★★★ is open-ended. The chapters used are in brackets.

1. ★ Diagonalising a symmetric matrix. Let S = \begin{pmatrix} 3 & 1 \\ 1 & 3 \end{pmatrix}. Compute its spectral decomposition, then S^{-1} and S^{10} without multiplying matrices (ch. 1).

   Solution: the eigenvalues are 4 and 2, with u_1 = (1, 1)/\sqrt2 and u_2 = (1, -1)/\sqrt2. Functions of S act on the eigenvalues: S^{-1} = \tfrac14 u_1 u_1^\top + \tfrac12 u_2 u_2^\top = \frac18 \begin{pmatrix} 3 & -1 \\ -1 & 3 \end{pmatrix}, \qquad S^{10} = 4^{10} u_1 u_1^\top + 2^{10} u_2 u_2^\top = \begin{pmatrix} 524800 & 523776 \\ 523776 & 524800 \end{pmatrix}.

2. ★ Three tests of positivity. Classify each matrix as positive definite, positive semidefinite or indefinite, with a different test each time (ch. 1): A = \begin{pmatrix} 2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 2 \end{pmatrix}, \qquad B = \begin{pmatrix} 1 & 2 \\ 2 & 3 \end{pmatrix}, \qquad C = \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}.

   Solution: the leading minors of A are 2, 3 and 4, so A is positive definite (Sylvester).

   The determinant of B is -1, so its eigenvalues have opposite signs; indeed x = (2, -1) gives x^\top B x = -1. Finally C = v v^\top with v = (1, 1) is a Gram matrix, hence PSD, with eigenvalues 2 and 0, so it is not definite.

3. ★ A regression through the origin. Fit y = a t to the points (1, 2), (2, 3), (3, 5) by least squares. Give \hat a, the residual, the projection matrix and the residual sum of squares, and check Theorem 2.4 (ch. 2).

   Solution: with t = (1, 2, 3) and y = (2, 3, 5), \hat a = t^\top y / t^\top t = 23/14 \approx 1.643.

   The residual y - \hat a t = (5, -4, 1)/14 is orthogonal to t, since 5 - 8 + 3 = 0. The projection matrix is t t^\top / 14, and the residual sum of squares is 42/196 = 3/14 \approx 0.214.

4. ★ Gradients of matrix functions. Compute (a) the gradient of \|x\| for x \neq 0; (b) the gradient of W \mapsto \frac12 \|XW - Y\|_F^2; (c) the gradient of f(K) = \log\det K - \operatorname{tr}(KS) over positive definite K, for a positive definite S, and its maximiser (ch. 4).

   Solution: (a) x / \|x\|, by the chain rule on \sqrt{x^\top x}. (b) X^\top (XW - Y), from Proposition 4.2 with B = I. (c) K^{-1} - S, which vanishes at K = S^{-1}.

   Along a line, t \mapsto \log\det(K + tH) has second derivative -\operatorname{tr}\big((K^{-1}H)^2\big) \le 0, so f is concave and S^{-1} is the maximiser.

   Interpretation: for centred observations x_i with S = \frac1n \sum_i x_i x_i^\top, the Gaussian log-likelihood of the precision matrix K is \frac n2 f(K) up to a constant. The maximum-likelihood covariance is therefore the empirical one, K^{-1} = S.

5. ★ Parseval and \pi^4/90. Compute the Fourier series of f(t) = |t| on (-\pi, \pi), extended with period 2\pi, and use Parseval to compute \sum_{n \ge 1} 1/n^4 (ch. 5).

   Solution: f is even, so b_n = 0, and a_0 = \pi.

   An integration by parts gives a_n = 2\big((-1)^n - 1\big)/(\pi n^2), that is -4/(\pi n^2) for odd n and 0 for even n \ge 2. Parseval (Example 5.7) gives \frac{2\pi^2}{3} = \frac{\pi^2}{2} + \sum_{n \text{ odd}} \frac{16}{\pi^2 n^4}, \qquad \text{so} \qquad \sum_{n \text{ odd}} \frac{1}{n^4} = \frac{\pi^4}{96} . The even terms are 1/16 of the full sum, so \sum_{n \ge 1} \frac{1}{n^4} = \frac{16}{15} \cdot \frac{\pi^4}{96} = \frac{\pi^4}{90} .

6. ★ Proximal operators by hand. Compute (a) the soft thresholding of v = (3, -0.5, 1.2) at level 1; (b) the projection of (1.5, -0.3) onto the square [0, 1]^2; (c) the prox of g = \frac12 \|x\|^2 with \gamma = 1 at (2, 4) (ch. 9).

   Solution: (a) (2, 0, 0.2); (b) (1, 0), by clipping each coordinate; (c) v / (1 + \gamma) = (1, 2).

7. ★★ Rank one and pseudo-inverse. For A = \begin{pmatrix} 1 & 1 \\ 1 & 1 \\ 1 & -1 \end{pmatrix}, compute the SVD, the best rank-one approximation and its error in both norms, and A^{+}; check that A^{+} A = I (ch. 3).

   Solution: A^\top A = \begin{pmatrix} 3 & 1 \\ 1 & 3 \end{pmatrix} has eigenvalues 4 and 2, so \sigma_1 = 2 and \sigma_2 = \sqrt2.

   The singular vectors are v_1 = (1, 1)/\sqrt2, v_2 = (1, -1)/\sqrt2, u_1 = (1, 1, 0)/\sqrt2 and u_2 = (0, 0, 1). Then A_1 = \begin{pmatrix} 1 & 1 \\ 1 & 1 \\ 0 & 0 \end{pmatrix}, \qquad A - A_1 = \begin{pmatrix} 0 & 0 \\ 0 & 0 \\ 1 & -1 \end{pmatrix}, \qquad A^{+} = \begin{pmatrix} 1/4 & 1/4 & 1/2 \\ 1/4 & 1/4 & -1/2 \end{pmatrix}. The error has operator and Frobenius norm \sqrt2 = \sigma_2, as Theorem 3.4 predicts, and a direct product gives A^{+} A = I.

8. ★★ Log-sum-exp. Let f(x) = \log \sum_i e^{x_i}. Show that \nabla f = p, the softmax vector, and that \nabla^2 f = \operatorname{diag}(p) - p p^\top. Deduce that f is convex and 1-smooth but not strongly convex (ch. 4, 6).

   Solution: \partial f / \partial x_i = e^{x_i} / \sum_j e^{x_j} = p_i, and \partial p_i / \partial x_j = p_i \delta_{ij} - p_i p_j.

   For every v, v^\top \nabla^2 f\, v = \sum_i p_i v_i^2 - (\sum_i p_i v_i)^2 is the variance of the values v_i under the weights p_i.

   It is non-negative, so f is convex, and at most \sum_i p_i v_i^2 \le \|v\|^2, so L \le 1.

   The Hessian kills the constant vector, because f(x + t\mathbf{1}) = f(x) + t is affine in t: no strong convexity. The softmax cross-entropy loss is f minus a linear function, with the same properties.

9. ★★ The best constant step on a quadratic. Let f(x) = \frac12 x^\top A x with \mu I \preceq A \preceq L I.

   Show that gradient descent contracts the error by \max(|1 - \eta\mu|, |1 - \eta L|) per step, and that \eta = 2/(L + \mu) minimises this factor. Compare with the step 1/L for \kappa = 100 and an accuracy of 10^{-3} (ch. 1, 7).

   Solution: x_{k+1} = (I - \eta A) x_k, and in the eigenbasis each component is multiplied by 1 - \eta \lambda_i. On [\mu, L], the maximum of |1 - \eta\lambda| is reached at an endpoint; it is smallest when 1 - \eta\mu = \eta L - 1, giving the factor (\kappa - 1)/(\kappa + 1).

   For \kappa = 100, the factor is 99/101 \approx 0.980, and \ln 1000 / |\ln 0.980| \approx 346 steps reach 10^{-3}. With \eta = 1/L the factor is 0.99 and about 688 steps are needed: the optimal step halves the work.

10. ★★ Projection onto a half-space. Project y onto \{x : a^\top x \le c\} with the KKT conditions, then compute it for a = (1, 2), c = 2 and y = (2, 2), with the multiplier and the dual value (ch. 8).

    Solution: stationarity of \frac12 \|x - y\|^2 + \lambda (a^\top x - c) gives x = y - \lambda a. If a^\top y \le c, then \lambda = 0 and x = y. Otherwise the constraint is active, and \lambda = \frac{a^\top y - c}{\|a\|^2}, \qquad x^* = y - \frac{\max(0,\ a^\top y - c)}{\|a\|^2}\, a . Here a^\top y = 6, so \lambda = 0.8 and x^* = (1.2, 0.4), on the boundary. The dual function g(\lambda) = \lambda (a^\top y - c) - \frac{\lambda^2}{2} \|a\|^2 is maximal at the same \lambda, with value 1.6 = \frac12 \|y - x^*\|^2.

11. ★★ Power iteration (Python). Build a symmetric 50 \times 50 matrix with \lambda_1 = 10 and \lambda_2 = 9, run power iteration and check the rate of Proposition 10.4 (ch. 1, 10).

    Solution:

    import numpy as np
    rng = np.random.default_rng(0)
    Q = np.linalg.qr(rng.standard_normal((50, 50)))[0]       # random orthonormal basis
    lam = np.r_[10.0, 9.0, rng.uniform(0, 5, 48)]            # lambda_2 / lambda_1 = 0.9
    S = Q @ np.diag(lam) @ Q.T
    u1 = Q[:, 0]
    x = rng.standard_normal(50)
    for k in range(1, 201):
        x = S @ x
        x /= np.linalg.norm(x)                               # power iteration
        if k in (50, 100, 200):
            sin = np.linalg.norm(x - (x @ u1) * u1)          # sine of the angle to u_1
            print(f"k = {k:3d}   sin = {sin:.1e}   ratio to 0.9^k = {sin / 0.9 ** k:.2f}")
    print(f"Rayleigh quotient {x @ S @ x:.6f}")

    Expected: the sine falls from 3.7 \times 10^{-3} at k = 50 to 5.1 \times 10^{-10} at k = 200, and its ratio to 0.9^k stays at 0.72. The Rayleigh quotient returns \lambda_1 = 10 to six decimals.

12. ★★ Gradient descent versus acceleration (Python). On a quadratic in dimension 100 with Hessian eigenvalues spread from 1 to 1000, count the iterations that gradient descent and Nesterov’s method need to divide f by 10^6 (ch. 7, 10).

    Solution:

    import numpy as np
    d = 100
    h = np.logspace(0, 3, d)               # Hessian eigenvalues from mu = 1 to L = 1000
    L, mu = h.max(), h.min()
    x0 = np.ones(d)
    f = lambda x: 0.5 * np.sum(h * x ** 2)  # minimum value 0 at x = 0
    
    def iterations(method, tol=1e-6):
        x, y = x0.copy(), x0.copy()
        beta = (np.sqrt(L) - np.sqrt(mu)) / (np.sqrt(L) + np.sqrt(mu))
        for k in range(1, 10 ** 6):
            if method == "gradient descent":
                x = x - h * x / L
            else:                            # Nesterov, constant momentum
                x_new = y - h * y / L
                y = x_new + beta * (x_new - x)
                x = x_new
            if f(x) <= tol * f(x0):
                return k
    
    for m in ("gradient descent", "Nesterov"):
        print(f"{m:16s} {iterations(m):6d} iterations")
    print(f"kappa = {L / mu:.0f}, sqrt(kappa) = {np.sqrt(L / mu):.1f}")

    Expected: 2685 iterations for gradient descent and 140 for Nesterov’s method, a ratio of 19. The theory predicts a ratio of order \sqrt\kappa \approx 32; the constants and the spread of the spectrum account for the difference.

13. ★★★ The dual of the lasso and a stopping rule. For P(x) = \frac12 \|Ax - b\|^2 + \lambda \|x\|_1, show that the dual problem is to maximise D(u) = -\frac12 \|u\|^2 - u^\top b subject to \|A^\top u\|_\infty \le \lambda. Deduce a computable bound on P(x) - P^* and use it to stop ISTA (ch. 2, 8, 9).

    Hints: write z = Ax - b as a constraint with multiplier u. The infimum over z of \frac12\|z\|^2 - u^\top z is -\frac12\|u\|^2.

    The infimum over x of \lambda\|x\|_1 + u^\top A x is 0 if \|A^\top u\|_\infty \le \lambda and -\infty otherwise, because \ell_\infty is the dual norm of \ell_1 (Definition 2.1).

    At the optimum u^* = Ax^* - b. Rescaling the current residual into the feasible set gives a dual point u and the bound P(x) - P^* \le P(x) - D(u).

    import numpy as np
    rng = np.random.default_rng(0)
    n, d, lam = 100, 200, 0.25                       # the lasso problem of ch. 9
    A = rng.standard_normal((n, d)) / np.sqrt(n)
    x_true = np.zeros(d)
    x_true[:5] = [3, -2, 2, -1.5, 1]
    b = A @ x_true + 0.05 * rng.standard_normal(n)
    L = np.linalg.norm(A, 2) ** 2
    
    def soft(v, t):
        return np.sign(v) * np.maximum(np.abs(v) - t, 0)
    
    def primal(x):
        return 0.5 * np.sum((A @ x - b) ** 2) + lam * np.abs(x).sum()
    
    def dual(u):                                      # valid when |A^T u|_inf <= lam
        return -0.5 * u @ u - u @ b
    
    x = np.zeros(d)
    for k in range(1, 2001):
        x = soft(x - A.T @ (A @ x - b) / L, lam / L)  # ISTA step
        r = A @ x - b
        u = r * min(1.0, lam / np.abs(A.T @ r).max())  # rescaled residual: dual feasible
        if k in (10, 100, 1000, 2000):
            print(f"k = {k:4d}   duality gap = {primal(x) - dual(u):.1e}")

    Expected: the gap is 1.2 after 10 iterations, 4.0 \times 10^{-6} after 100 and reaches rounding level, about 10^{-14}, by 1000. A stopping rule “gap below 10^{-6}” is a certificate, unlike a rule on the change of the iterates.

14. ★★★ The edge of stability. Train a small multilayer perceptron (two hidden layers of width 32, tanh) by full-batch gradient descent on half the mean squared error. Use 500 Gaussian inputs in \mathbb{R}^5 and targets \sin(w^\top x) + 0.1\,\varepsilon, with random w and standard Gaussian noise \varepsilon.

    Run 5000 steps with \eta \in \{0.1, 0.2, 0.4\}. Track the loss and the largest Hessian eigenvalue, the sharpness, and compare it with 2/\eta (ch. 6, 7, 10, 11).

    Hints and expected result: compute Hessian–vector products with torch.autograd.grad applied twice, and the sharpness with 20 steps of power iteration (Proposition 10.4) every 50 iterations. Each run takes a few seconds on a CPU.

    The sharpness starts between 2 and 3, rises, then hovers at 2/\eta (20, 10 and 5), as Cohen and coauthors (2021) report; for \eta = 0.1 the rise can take most of the 5000 steps. Meanwhile the loss keeps decreasing over the run, but not monotonically. With \eta = 0.01, the sharpness is still near 6 after 3000 steps, far from 2/\eta = 200: the regime needs large steps or long runs.

    Explain why Lemma 6.4 says nothing in this regime, and discuss what replaces it in the paper.

Projects and assessment

The refresher fortnight announces no assessment: per the master-mva.com planning, the refresher courses are optional and without registration, and mandatory for the I&I track, with no grade mentioned. This section therefore replaces the usual projects and written exam by a diagnostic self-check, to take before and after the four steps, and two small practice projects that prepare the paper studies common in MVA courses. None of it is official material.

Diagnostic self-check

Twelve short questions, about two minutes each, one or two per chapter. Cover the answer while you think.

1. Q1 — eigenvalues (ch. 1). What are the eigenvalues and an orthonormal eigenbasis of \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}?

   Answer: 1 and -1, with eigenvectors (1, 1)/\sqrt2 and (1, -1)/\sqrt2.

2. Q2 — positivity (ch. 1). Is \begin{pmatrix} 1 & 2 \\ 2 & 1 \end{pmatrix} positive semidefinite?

   Answer: no. Its determinant is -3 and its eigenvalues are 3 and -1; the vector (1, -1) gives the quadratic form -2.

3. Q3 — projection (ch. 2). Project (1, 2, 3) onto the line spanned by (1, 1, 1).

   Answer: (2, 2, 2); the residual (-1, 0, 1) is orthogonal to (1, 1, 1).

4. Q4 — norms (ch. 2). Give the operator norm, the Frobenius norm and the condition number of \operatorname{diag}(3, -4).

   Answer: 4, 5 and 4/3.

5. Q5 — gradients (ch. 0). Give the gradient and the Hessian of \|Ax - b\|^2, without the factor \frac12.

   Answer: 2 A^\top (Ax - b) and 2 A^\top A.

6. Q6 — SVD (ch. 3). What are the rank, the singular values and the nuclear norm of u v^\top with u = (1, 2, 2) and v = (3, 4)?

   Answer: rank 1, a single non-zero singular value \|u\| \|v\| = 3 \times 5 = 15; the nuclear, operator and Frobenius norms all equal 15.

7. Q7 — matrix calculus (ch. 4). What is the gradient of X \mapsto \operatorname{tr}(AXB)?

   Answer: A^\top B^\top, because the cyclic property of the trace gives \operatorname{tr}(AHB) = \operatorname{tr}(BAH) = \langle (BA)^\top, H \rangle_F .

8. Q8 — Fourier (ch. 5). What are the Fourier transforms of f(x - a) and of e^{2i\pi b x} f(x)?

   Answer: e^{-2i\pi a \xi} \hat f(\xi) and \hat f(\xi - b).

9. Q9 — convexity (ch. 6). Is f(x) = \max(x_1, x_2) + x_1^2 convex? Strongly convex?

   Answer: convex, as the sum of a maximum of linear functions and a convex quadratic. Not strongly convex: f(0, t) = \max(0, t) is piecewise linear in t.

10. Q10 — gradient descent (ch. 7). For f(x) = \frac12 (x_1^2 + 10 x_2^2), what is L, which constant steps converge, and how many steps of size 1/L divide f by 10^6?

    Answer: L = 10; every step \eta < 2/L = 0.2 converges. With \eta = 1/L, the first coordinate is multiplied by 0.9 per step and f by 0.81, so about 66 steps.

11. Q11 — duality (ch. 8). Write the dual of “minimise x over x \in \mathbb{R} subject to x \ge 1”.

    Answer: the Lagrangian x + \lambda(1 - x) has infimum \lambda if \lambda = 1 and -\infty otherwise. So d^* = 1 = p^*, with \lambda^* = 1 and an active constraint.

12. Q12 — prox (ch. 9). Compute the prox of 0.5\, |\cdot| at -0.3 and at 2.

    Answer: 0, since |-0.3| \le 0.5, and 1.5.

Practice projects

Two small projects, sized for the refresher fortnight. Each reproduces a figure of a real paper, the format of most MVA course projects.

1. [Solo · ≈ 6 h] Randomized SVD: a paper study. Reproduce, at a small scale, the main experiment of Halko, Martinsson and Tropp (2011).
   • Implement: the randomized range finder (a Gaussian test matrix with oversampling p, then QR) with q power iterations, followed by the SVD of the small projected matrix.
   • Data: the digits of scikit-learn (1797 \times 64) and a synthetic 2000 \times 1000 matrix with prescribed singular values, decaying like 1/j or like 0.9^j.
   • Success criterion: a plot of the ratio of the spectral error to \sigma_{k+1} against q = 0, 1, 2, 3, for p = 5 and p = 10 and the rank k = 20 on the synthetic matrix. The ratio lies between 1.2 and 2 without power iterations and is close to 1 from two of them on. Add a timing comparison with numpy.linalg.svd.
   • Deliverables: a notebook and a two-page note relating the curves to Theorem 3.4 and Proposition 10.4. Chapters 1, 3, 10.
2. [Group of 2 · ≈ 8 h each] First-order methods for the lasso: a paper study. Beck and Teboulle (2009) compare ISTA and FISTA on wavelet-regularised image deblurring. Reproduce their comparison on a one-dimensional deconvolution, then contrast it with sparse recovery.
   • Implement: ISTA, FISTA (Algorithm 10.7 with a prox step) and a reference solver (scikit-learn’s Lasso), all stopped by the duality gap of Exercise 13.
   • Data: 15 spikes on a signal of length 500, blurred by a Gaussian kernel of standard deviation 4 samples, with noise of standard deviation 0.01 and \lambda = 0.01 \|A^\top b\|_\infty. Then a Gaussian 200 \times 500 matrix with a 10-sparse signal, noise of standard deviation 0.05 and \lambda = 0.1 \|A^\top b\|_\infty.
   • Success criterion: on the deconvolution, a log–log plot of F(x_k) - F^* for k from 100 to 10 000, below the bounds of Algorithms 9.7 and 10.7. The slope of ISTA stays near -1; that of FISTA is -2 or steeper up to k = 1000, before a final linear phase.
   • Second criterion: on the Gaussian matrix, both methods reach 10^{-12} within a few hundred steps. Explain this linear rate by the strong convexity of the problem restricted to the support, and recover the support along a grid of \lambda.
   • Work split: one member writes the solvers and the duality gap, the other the experiments and the regularisation path; both write the report.
   • Deliverables: code, a four-page report and a ten-minute talk. Chapters 6–10.

Self-assessment grid

References

The refresher fortnight publishes no programme and no reference list (master-mva.com planning, 2026–2027). The works below are standard references at the level of this page, followed by the papers cited in it.

Linear algebra and matrix calculus

• G. Strang, Linear Algebra and Learning from Data, Wellesley-Cambridge Press, 2019.
• R. A. Horn, C. R. Johnson, Matrix Analysis, 2nd ed., Cambridge University Press, 2013.
• L. N. Trefethen, D. Bau III, Numerical Linear Algebra, SIAM, 1997.
• G. H. Golub, C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, 2013.
• J. R. Magnus, H. Neudecker, Matrix Differential Calculus with Applications in Statistics and Econometrics, 3rd ed., Wiley, 2019.
• K. B. Petersen, M. S. Pedersen, The Matrix Cookbook, Technical University of Denmark, 2012.
• C. Eckart, G. Young, “The approximation of one matrix by another of lower rank”, Psychometrika 1(3), 1936.
• L. Mirsky, “Symmetric gauge functions and unitarily invariant norms”, Quarterly Journal of Mathematics 11(1), 1960.

Analysis and Fourier analysis

• H. Brezis, Functional Analysis, Sobolev Spaces and Partial Differential Equations, Springer, 2011.
• E. M. Stein, R. Shakarchi, Fourier Analysis: An Introduction, Princeton University Press, 2003.
• S. Mallat, A Wavelet Tour of Signal Processing: The Sparse Way, 3rd ed., Academic Press, 2009.
• J. W. Cooley, J. W. Tukey, “An algorithm for the machine calculation of complex Fourier series”, Mathematics of Computation 19(90), 1965.

Optimisation

• S. Boyd, L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004.
• Y. Nesterov, Lectures on Convex Optimization, 2nd ed., Springer, 2018.
• J. Nocedal, S. J. Wright, Numerical Optimization, 2nd ed., Springer, 2006.
• A. Beck, First-Order Methods in Optimization, SIAM, 2017.
• H. H. Bauschke, P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, 2nd ed., Springer, 2017.
• S. Bubeck, “Convex optimization: algorithms and complexity”, Foundations and Trends in Machine Learning 8(3–4), 2015.
• L. Bottou, F. E. Curtis, J. Nocedal, “Optimization methods for large-scale machine learning”, SIAM Review 60(2), 2018.
• N. Parikh, S. Boyd, “Proximal algorithms”, Foundations and Trends in Optimization 1(3), 2014.
• H. W. Kuhn, A. W. Tucker, “Nonlinear programming”, Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability, University of California Press, 1951.
• H. Robbins, S. Monro, “A stochastic approximation method”, Annals of Mathematical Statistics 22(3), 1951.
• A. Nemirovski, A. Juditsky, G. Lan, A. Shapiro, “Robust stochastic approximation approach to stochastic programming”, SIAM Journal on Optimization 19(4), 2009.
• Y. Nesterov, “A method for solving the convex programming problem with convergence rate O(1/k²)”, Soviet Mathematics Doklady 27, 1983.
• A. Beck, M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems”, SIAM Journal on Imaging Sciences 2(1), 2009.

Research directions (ch. 11)

• N. Halko, P.-G. Martinsson, J. A. Tropp, “Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions”, SIAM Review 53(2), 2011.
• P.-G. Martinsson, J. A. Tropp, “Randomized numerical linear algebra: foundations and algorithms”, Acta Numerica 29, 2020.
• E. J. Candès, B. Recht, “Exact matrix completion via convex optimization”, Foundations of Computational Mathematics 9(6), 2009.
• S. Burer, R. D. C. Monteiro, “A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization”, Mathematical Programming 95(2), 2003.
• R. Ge, J. D. Lee, T. Ma, “Matrix completion has no spurious local minimum”, NeurIPS 2016.
• S. Gunasekar, B. Woodworth, S. Bhojanapalli, B. Neyshabur, N. Srebro, “Implicit regularization in matrix factorization”, NeurIPS 2017.
• S. Arora, N. Cohen, W. Hu, Y. Luo, “Implicit regularization in deep matrix factorization”, NeurIPS 2019.
• E. J. Hu, Y. Shen, P. Wallis, Z. Allen-Zhu, Y. Li, S. Wang, L. Wang, W. Chen, “LoRA: low-rank adaptation of large language models”, ICLR 2022.
• D. Soudry, E. Hoffer, M. S. Nacson, S. Gunasekar, N. Srebro, “The implicit bias of gradient descent on separable data”, Journal of Machine Learning Research 19(70), 2018.
• J. Cohen, S. Kaur, Y. Li, J. Z. Kolter, A. Talwalkar, “Gradient descent on neural networks typically occurs at the edge of stability”, ICLR 2021.
• B. T. Polyak, “Gradient methods for the minimisation of functionals”, USSR Computational Mathematics and Mathematical Physics 3(4), 1963.
• H. Karimi, J. Nutini, M. Schmidt, “Linear convergence of gradient and proximal-gradient methods under the Polyak–Łojasiewicz condition”, ECML PKDD 2016.
• Y. Drori, M. Teboulle, “Performance of first-order methods for smooth convex minimization: a novel approach”, Mathematical Programming 145, 2014.
• A. B. Taylor, J. M. Hendrickx, F. Glineur, “Smooth strongly convex interpolation and exact worst-case performance of first-order methods”, Mathematical Programming 161, 2017.
• M. Schmidt, N. Le Roux, F. Bach, “Minimizing finite sums with the stochastic average gradient”, Mathematical Programming 162, 2017.
• R. Johnson, T. Zhang, “Accelerating stochastic gradient descent using predictive variance reduction”, NeurIPS 2013.
• A. Defazio, F. Bach, S. Lacoste-Julien, “SAGA: a fast incremental gradient method with support for non-strongly convex composite objectives”, NeurIPS 2014.
• F. Bach, E. Moulines, “Non-strongly-convex smooth stochastic approximation with convergence rate O(1/n)”, NeurIPS 2013.
• D. P. Kingma, J. Ba, “Adam: a method for stochastic optimization”, ICLR 2015.
• V. Gupta, T. Koren, Y. Singer, “Shampoo: preconditioned stochastic tensor optimization”, ICML 2018.
• S. V. Venkatakrishnan, C. A. Bouman, B. Wohlberg, “Plug-and-play priors for model based reconstruction”, IEEE Global Conference on Signal and Information Processing (GlobalSIP), 2013.
• E. K. Ryu, J. Liu, S. Wang, X. Chen, Z. Wang, W. Yin, “Plug-and-play methods provably converge with properly trained denoisers”, ICML 2019.
• S. Hurault, A. Leclaire, N. Papadakis, “Gradient step denoiser for convergent plug-and-play”, ICLR 2022.