Foundations of Distributed and Large-Scale Computing Optimization

Contents

• Introduction and overview
• 0. Review: convex analysis and gradient descent
• 1. Large-scale optimization: models, costs and accelerated gradients
• 2. Proximity operators
• 3. Proximal gradient algorithms: forward–backward and FISTA
• 4. Splitting algorithms: Douglas–Rachford, ADMM and primal–dual methods
• 5. Majorization–Minimization methods
• 6. Stochastic and block-coordinate methods for large data sets
• 7. Parallel optimization
• 8. Distributed optimization strategies
• 9. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student with a good M1 in applied mathematics, computer science or signal processing who can minimize a smooth convex function by gradient descent and now faces problems with millions of unknowns, non-smooth regularizers, or data spread over several machines.

Empirical first: on a denoising problem with 2000 unknowns and an edge-preserving penalty, plain gradient descent still has an objective error of 25.6 after 50 iterations. A Majorization–Minimization method of about the same cost per iteration is at 2.8 \times 10^{-7} (ch. 5). On a ring of 10 machines holding different data, decentralized gradient descent stalls at distance 0.42 from the solution forever, while gradient tracking reaches machine precision (ch. 8).

This course explains both observations. It builds the convex-analytic toolbox (proximity operators, monotone operators, majorants) behind the algorithms that solve very large problems. It then shows how the same algorithms split across the cores of one machine or across a network of machines.

The course shares tools with Convex Optimization and Applications in Machine Learning, taught in the same semester, and prepares the imaging and large-model courses of semester 2, such as Deep Learning for Image Restoration and Synthesis and Training and Deploying Large-Scale Models.

What you will be able to do

• Compute the proximity operator of the usual regularizers and constraints, and combine them with the rules of proximal calculus.
• Prove the convergence rates of gradient descent, forward–backward and FISTA from the descent lemma.
• Choose between forward–backward, Douglas–Rachford, ADMM and primal–dual methods from the structure of a problem.
• Build a tangent majorant for a given cost (Lipschitz, half-quadratic, Jensen) and derive the Majorization–Minimization algorithm it defines.
• Implement the subspace Majorization–Minimization method 3MG and explain why it scales to large images.
• Analyze stochastic gradient, variance-reduced and randomized coordinate methods on finite sums.
• Reformulate a sum of functions in a product space and run the parallel proximal algorithm PPXA.
• Predict the speedup of a parallel method with Amdahl’s law, the cost of an all-reduce and the damping of parallel coordinate steps.
• Design a distributed method (consensus ADMM, decentralized gradient, gradient tracking) for data held by several agents, and explain the bias of naive schemes.
• Implement these methods in NumPy and PyTorch and check them against a reference solution.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com states as objective to introduce the theoretical background that makes it possible to develop efficient algorithms for the large-scale problems the course addresses, taking advantage of modern multicore or distributed computing architectures.

The description adds that the course is mainly focused on nonlinear optimization tools for convex problems. It presents proximal tools, splitting techniques and Majorization–Minimization strategies, which the description calls very popular for processing massive data sets, and it illustrates these methods on various application examples.

The listed instructor is Emilie Chouzenoux (CentraleSupélec & INRIA Saclay). The organisation announced is 7 sessions of 3 h plus a 2 h exam, and the assessment is two lab reports plus an exam. The course page also links a presentation file; these notes are an independent reconstruction and do not reproduce it. The Université Paris-Saclay course detail gives only the type of assessment, with no programme of its own.

Every official topic is taught on this page. The extensions beyond the syllabus are the review of ch. 0, the lower bound and acceleration theory of ch. 1, the stochastic and variance-reduced methods of ch. 6, federated learning in ch. 8, and the self-study research chapter 9.

Where this course fits in the programme

Official prerequisites: not displayed. In practice the course uses M1 linear algebra (eigenvalues, symmetric matrices), differential calculus and the basics of convex functions; ch. 0 reviews what is needed. Chapter 6 also uses expectations and variances.

• Refreshers — Linear Algebra, Analysis and Optimization for ch. 0, Probability and Statistics for the stochastic methods of ch. 6, and Scientific Python, NumPy and PyTorch for the code.
• Related courses of semester 1 — Convex Optimization and Applications in Machine Learning (duality and interior-point methods, complementary to the first-order methods here), Computational Statistics (EM and stochastic approximation, met as MM and SGD), Introduction to Probabilistic Graphical Models (EM), Deep Learning (stochastic training), Optimal Transport (entropic solvers as alternating projections) and Image Denoising (variational denoisers).
• Semester 2 courses that build on it — Deep Learning for Image Restoration and Synthesis (variational and hybrid methods for inverse problems), Inverse Problems and Imaging (regularized least squares for ill-posed problems), Training and Deploying Large-Scale Models (data parallelism and all-reduce), Algorithms for Network Optimization and Management (optimization problems of communication networks), Graphs in Machine Learning (graph Laplacians, the matrices of gossip), Sequential Learning (online convex optimization) and Kernel Methods (large-scale kernel solvers).

Course outline

• 0. Review: convex analysis and gradient descent — subdifferentials, smoothness and strong convexity, the rates of gradient descent, conjugates and Fenchel duality.
• 1. Large-scale optimization — the composite model, what second-order methods cost, the lower bound of first-order methods, acceleration, estimating Lipschitz constants and preconditioning.
• 2. Proximity operators — definition, a catalogue, firm nonexpansiveness, Moreau decomposition and envelope, proximal calculus.
• 3. Proximal gradient algorithms — forward–backward, its O(1/k) rate, FISTA and its O(1/k^2) rate, variable metrics.
• 4. Splitting algorithms — averaged operators, Douglas–Rachford, ADMM, Chambolle–Pock and Condat–Vũ.
• 5. Majorization–Minimization — tangent majorants, half-quadratic and EM, convergence, the subspace method 3MG.
• 6. Stochastic and block-coordinate methods — stochastic gradient, mini-batches, SVRG and SAGA, randomized coordinate descent.
• 7. Parallel optimization — speedup and communication models, PPXA in a product space, parallel coordinate descent, data parallelism and asynchrony.
• 8. Distributed optimization strategies — consensus ADMM, gossip, decentralized gradient, gradient tracking, federated averaging.
• 9. Research directions and open problems — a self-study tour of current work.

Session plan

24 h in class = 8 sessions of 3 h. The volume of 24 h of lectures comes from the Université Paris-Saclay programme table. master-mva.com organises the course as 7 sessions of 3 h plus a 2 h exam. The split below is a proposal: it keeps seven 3-hour lectures, uses the remaining 3 h of the official volume for an eighth session, and leaves the written exam outside the plan.

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

Self-study, outside the 24 h (about 4 h in total): ch. 9 Research directions and open problems (≈ 4 h).

How to study this module

• Read ch. 0 before session 1 if subdifferentials or conjugate functions feel distant: every later proof uses them.
• Follow one thread through the course: almost every algorithm is a fixed-point iteration of an averaged operator (ch. 4), and almost every proof is the descent lemma of ch. 0 plus one inequality.
• After each session, do the exercises listed in its row, then run the code of the chapter and change one parameter (step size, number of workers, graph) to see the effect predicted by the theory.
• The code uses NumPy, SciPy, scikit-learn, NetworkX and PyTorch (CPU) on synthetic data; every block runs in a few seconds on a laptop.
• The solo projects of the last section are sized like the two official lab reports; start the first one after session 3.

0. Review: convex analysis and gradient descent

Every algorithm of the course minimizes a convex function that may be non-smooth or take the value +\infty. This chapter collects the M1-level facts the proofs rely on: subgradients, the descent lemma, the rates of gradient descent and Fenchel duality.

The Linear Algebra, Analysis and Optimization refresher covers the same ground at a slower pace. Throughout, the space is \mathbb{R}^N with the Euclidean inner product \langle \cdot, \cdot \rangle; all results hold in Hilbert spaces (Bauschke and Combettes, 2011).

Convex functions and subdifferentials

Indicators turn constraints into terms of the objective: minimizing f over C is minimizing f + \iota_C over the whole space.

Smoothness, strong convexity and gradient descent

The descent lemma says that a smooth function lies below a parabola tangent at any point. Minimizing that parabola is a gradient step: it is the first Majorization–Minimization algorithm of the course (ch. 5).

To reach a relative accuracy \varepsilon on \|x_k - x^\star\|^2, the linear rate needs about \kappa \ln(1/\varepsilon) iterations. Chapter 1 shows that this dependence on \kappa is not optimal.

Conjugates and duality

The dual variable u lives in the space of Kx. Primal–dual algorithms (ch. 4) update x and u together and only ever multiply by K and K^\top.

1. Large-scale optimization: models, costs and accelerated gradients

What makes an optimization problem “large”, and which algorithms survive at that scale? This chapter fixes the model problem of the course, shows why second-order methods are out of reach, and establishes what first-order methods can and cannot achieve.

The model problem

Most problems of the course have the composite form \min_{x \in \mathbb{R}^N} \; F(x) = \sum_{j=1}^{J} f_j(L_j x), with each f_j \in \Gamma_0 simple in some sense (smooth, or with an explicit proximity operator) and each L_j a linear operator that can be applied quickly.

“Large” means that N, n or both are too large for any method whose cost per iteration grows faster than linearly with the size of the data.

What first-order methods can achieve

Gradient descent reaches O(1/k) (Theorem 0.7), while the bound allows O(1/k^2). Nesterov (1983) closed the gap with a momentum term.

The gain is large on ill-conditioned problems. With \kappa = 10^4 and \varepsilon = 10^{-6}, the two iteration counts are \kappa \ln(1/\varepsilon) \approx 1.4 \times 10^5 \ \text{ (gradient descent)}, \qquad \sqrt{\kappa} \ln(1/\varepsilon) \approx 1.4 \times 10^3 \ \text{ (accelerated)}.

Estimating the Lipschitz constant

Every method above needs L, and a step larger than 2/L can diverge. For f(x) = \frac12 \|Ax - b\|^2, the constant is L = \|A\|^2, the largest eigenvalue of A^\top A, which is too costly to compute exactly at scale.

Preconditioning

Ill-conditioning slows every first-order method. A preconditioner changes the geometry in which the gradient step is taken.

A diagonal preconditioner costs nothing to invert, and its step updates all coordinates independently. This is the first example of a separable majorant, which makes an algorithm parallel (ch. 5 and ch. 7).

The code below estimates L by power iteration on an ill-conditioned least-squares problem (singular values of the design spread over two decades), then compares gradient descent and Algorithm 1.5.

import numpy as np

rng = np.random.default_rng(0)
n, d = 300, 200
A = rng.standard_normal((n, d)) @ np.diag(np.logspace(0, -2, d))  # ill-conditioned
x_true = rng.standard_normal(d)
b = A @ x_true + 0.01 * rng.standard_normal(n)

# Lipschitz constant of grad f, f(x) = 0.5 ||Ax - b||^2, by power iteration
v = rng.standard_normal(d)
for _ in range(100):
    v = A.T @ (A @ v)
    v /= np.linalg.norm(v)
L = v @ (A.T @ (A @ v))
print(f"power iteration L = {L:.4f}, exact ||A||^2 = {np.linalg.norm(A, 2)**2:.4f}")

f = lambda x: 0.5 * np.sum((A @ x - b) ** 2)
f_star = f(np.linalg.lstsq(A, b, rcond=None)[0])
grad = lambda x: A.T @ (A @ x - b)

x = np.zeros(d)                           # gradient descent, step 1/L
y, z, t = np.zeros(d), np.zeros(d), 1.0   # Nesterov (FISTA with g = 0)
for k in range(1, 501):
    x = x - grad(x) / L
    z_new = y - grad(y) / L
    t_new = (1 + np.sqrt(1 + 4 * t * t)) / 2
    y = z_new + (t - 1) / t_new * (z_new - z)
    z, t = z_new, t_new
    if k in (10, 100, 500):
        print(f"k={k:4d}  GD gap {f(x) - f_star:.3e}   Nesterov gap {f(z) - f_star:.3e}")

Power iteration recovers L = 348.1558 to all printed digits. After 500 iterations the objective gaps are \text{gradient descent: } 2.81, \qquad \text{accelerated: } 7.6 \times 10^{-3} . At k = 100 the gaps are already 15.3 and 0.86.

2. Proximity operators

Gradient steps need differentiability, and the \ell_1 norm, total variation or a constraint are not differentiable. The proximity operator replaces the gradient step on such a term by a small, often explicit, minimization problem.

Definition and catalogue

The inclusion reads p = (I + \gamma \partial f)^{-1}(v): the prox is an implicit (backward) gradient step, where the subgradient is taken at the arrival point.

The proximal point algorithm x_{k+1} = \mathrm{prox}_{\gamma f}(x_k) converges to a minimizer of f whenever one exists (Rockafellar, 1976). It is rarely practical on its own, since each step is as hard as the original problem, but it is the template of every algorithm of chapters 3 and 4.

Firm nonexpansiveness and Moreau’s identities

When f is a norm, f^* is the indicator of the dual unit ball (Example 0.9), and scaling an indicator changes nothing. The prox of a norm is therefore the identity minus a projection: soft thresholding is v - P_{[-\gamma, \gamma]^N}(v). This identity also gives the prox of a conjugate, which the primal–dual methods of ch. 4 need.

Proximal calculus

Rule 3 covers \ell_1 regularization in an orthonormal wavelet or cosine basis. It does not cover total variation, because the discrete gradient D is not semi-orthogonal. In general, the prox of a sum is not the composition of the proxes: that is why splitting algorithms exist (ch. 4).

The code checks soft thresholding, the Moreau decomposition for the \ell_1 norm and firm nonexpansiveness on random pairs.

import numpy as np

def prox_l1(v, gamma):
    """prox of gamma * ||.||_1: soft thresholding"""
    return np.sign(v) * np.maximum(np.abs(v) - gamma, 0.0) + 0.0

def proj_linf(v, r):
    """projection onto the l-infinity ball of radius r"""
    return np.clip(v, -r, r)

rng = np.random.default_rng(1)
v = np.array([3.0, -0.4, 1.2, -2.5, 0.05])
print("prox_l1(v, 1)   =", prox_l1(v, 1.0))
# Moreau decomposition: v = prox_{gamma f}(v) + gamma * prox_{f*/gamma}(v/gamma),
# with f = ||.||_1, f* = indicator of the unit l-infinity ball
print("Moreau residual =", np.abs(v - prox_l1(v, 1.0) - proj_linf(v, 1.0)).max())

# firm nonexpansiveness: ||p - q||^2 <= <p - q, u - w>
worst = np.inf
for _ in range(10000):
    u, w = rng.standard_normal(5) * 3, rng.standard_normal(5) * 3
    p, q = prox_l1(u, 0.7), prox_l1(w, 0.7)
    worst = min(worst, (p - q) @ (u - w) - (p - q) @ (p - q))
print(f"min of <p-q,u-w> - ||p-q||^2 over 10^4 pairs: {worst:.2e}")

The output is (2, 0, 0.2, -1.5, 0), a Moreau residual of exactly 0, and a minimum of -7 \times 10^{-15} for the firm-nonexpansiveness gap, which is zero up to rounding.

3. Proximal gradient algorithms: forward–backward and FISTA

Many problems split into a smooth data term and a non-smooth but proximable regularizer. This chapter combines a gradient step on the first with a proximal step on the second, and proves the rates of the plain and accelerated versions.

Throughout the chapter, F = f + g with f convex and L-smooth, g \in \Gamma_0(\mathbb{R}^N), and a minimizer x^\star exists.

Forward–backward splitting

Special cases: g = 0 gives gradient descent; f = 0 gives the proximal point algorithm; g = \iota_C gives projected gradient; g = \lambda \|\cdot\|_1 gives ISTA, the iterative shrinkage-thresholding algorithm (Daubechies, Defrise and De Mol, 2004).

The norm of x_k - x_{k+1}, divided by \gamma_k, therefore measures optimality and gives a stopping criterion.

Acceleration: FISTA

FISTA is not monotone, and the convergence of its iterates is subtle; Chambolle and Dossal (2015) prove it for the slightly modified sequence t_k = (k + a - 1)/a with a > 2.

Variable metric and line search

When L is unknown, backtracking multiplies a trial L by a factor \eta > 1 until the descent inequality of Lemma 0.6 holds between y_k and the new point. The rates are kept, with L replaced by \eta L in the worst case (Beck and Teboulle, 2009).

import numpy as np
from sklearn.linear_model import Lasso

rng = np.random.default_rng(0)
n, d, s = 200, 500, 10
A = rng.standard_normal((n, d)) / np.sqrt(n)
x_true = np.zeros(d)
x_true[rng.choice(d, s, replace=False)] = rng.choice([-1.0, 1.0], s) * rng.uniform(1, 2, s)
b = A @ x_true + 0.01 * rng.standard_normal(n)
lam = 0.02

F = lambda x: 0.5 * np.sum((A @ x - b) ** 2) + lam * np.abs(x).sum()
soft = lambda v, g: np.sign(v) * np.maximum(np.abs(v) - g, 0.0)
L = np.linalg.norm(A, 2) ** 2

# reference solution: scikit-learn minimises (1/2n)||b - Ax||^2 + alpha ||x||_1
ref = Lasso(alpha=lam / n, fit_intercept=False, tol=1e-12, max_iter=100000).fit(A, b).coef_
F_star = F(ref)

x = np.zeros(d)                          # ISTA = forward-backward
z, y, t = np.zeros(d), np.zeros(d), 1.0  # FISTA
for k in range(1, 301):
    x = soft(x - A.T @ (A @ x - b) / L, lam / L)
    z_new = soft(y - A.T @ (A @ y - b) / L, lam / L)
    t_new = (1 + np.sqrt(1 + 4 * t * t)) / 2
    y = z_new + (t - 1) / t_new * (z_new - z)
    z, t = z_new, t_new
    if k in (10, 50, 300):
        print(f"k={k:5d}  ISTA {F(x) - F_star:.2e}   FISTA {F(z) - F_star:.2e}")

supp = lambda v: set(np.flatnonzero(np.abs(v) > 1e-6))
print("true support recovered by FISTA:", supp(x_true) <= supp(z), "| nonzeros:", len(supp(z)))
print(f"relative error ||x - x_true|| / ||x_true|| = {np.linalg.norm(z - x_true) / np.linalg.norm(x_true):.3f}")

The objective gaps are:

The final estimate has 38 non-zeros: the 10 true ones and 28 spurious ones, all below 0.02 in absolute value. Its relative error is 0.024. Both methods end up converging linearly, because the lasso is strongly convex on the support of the solution; FISTA reaches that regime much earlier.

4. Splitting algorithms: Douglas–Rachford, ADMM and primal–dual methods

Forward–backward needs one smooth term. Total variation \lambda\|Dx\|_1, a sum of two constraints, or a non-smooth data term break this assumption. This chapter builds algorithms that only use the proxes of the individual terms and products by linear operators, from one fixed-point theorem.

Averaged operators

By Proposition 4.2, the forward–backward operator \mathrm{prox}_{\gamma g} \circ (I - \gamma \nabla f) is averaged for \gamma \in (0, 2/L), and its fixed points are the minimizers (Proposition 3.2). Theorem 4.3 thus proves the convergence of the iterates claimed in Theorem 3.4.

Douglas–Rachford splitting

Douglas–Rachford uses no gradient and no Lipschitz constant: any \gamma > 0 works, although \gamma strongly affects the speed (Example 7.7).

ADMM

The alternating direction method of multipliers handles a linear coupling between two blocks of variables: \min_{x, z} \; f(x) + g(z) \quad \text{subject to} \quad A x + B z = c .

Boyd et al. (2011) recommend stopping when both the primal residual \|A x_k + B z_k - c\| and the dual residual \rho \|A^\top B (z_{k} - z_{k-1})\| are small.

Primal–dual methods

The problem \min_x g(x) + h(Kx) is equivalent to the saddle-point problem \min_x \max_u \; g(x) + \langle K x, u \rangle - h^*(u), by Definition 0.8 applied to h = h^{**}. Primal–dual methods alternate a prox step in x and one in u, and never invert K.

Condat (2013) and Vũ (2013) add a smooth term f with an L-Lipschitz gradient to the primal step: x_{k+1} = \mathrm{prox}_{\tau g}\big(x_k - \tau (\nabla f(x_k) + K^\top u_{k+1})\big), \qquad \text{with} \quad \frac{1}{\tau} - \sigma \|K\|^2 > \frac{L}{2} . The result handles f + g + h \circ K with only explicit operations: the most flexible of the splitting algorithms.

import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import splu

rng = np.random.default_rng(0)
N, lam = 500, 1.0
x_clean = np.repeat([0.0, 2.0, -1.0, 1.0, 0.0], N // 5)       # piecewise constant
y = x_clean + 0.5 * rng.standard_normal(N)
D = sp.diags([-np.ones(N - 1), np.ones(N - 1)], [0, 1], shape=(N - 1, N)).tocsr()
F = lambda x: 0.5 * np.sum((x - y) ** 2) + lam * np.abs(D @ x).sum()

def chambolle_pock(iters, tau=0.49, sigma=0.49):        # tau * sigma * ||D||^2 < 1
    x, xbar, u = y.copy(), y.copy(), np.zeros(N - 1)
    for _ in range(iters):
        u = np.clip(u + sigma * (D @ xbar), -lam, lam)    # prox of sigma h*, h = lam ||.||_1
        x_new = (x - tau * (D.T @ u) + tau * y) / (1 + tau)  # prox of tau g, g = 0.5||. - y||^2
        xbar, x = 2 * x_new - x, x_new
    return x

def admm(iters, rho=1.0):
    solve = splu((sp.eye(N) + rho * D.T @ D).tocsc()).solve   # factorised once
    z, w = np.zeros(N - 1), np.zeros(N - 1)
    for _ in range(iters):
        x = solve(y + rho * (D.T @ (z - w)))
        Dx = D @ x
        z = np.sign(Dx + w) * np.maximum(np.abs(Dx + w) - lam / rho, 0.0)
        w = w + Dx - z
    return x

F_star = F(chambolle_pock(20000))
for k in (50, 200, 1000):
    print(f"k={k:5d}  Chambolle-Pock {F(chambolle_pock(k)) - F_star:.2e}   ADMM {F(admm(k)) - F_star:.2e}")
x = chambolle_pock(20000)
print(f"RMSE noisy {np.sqrt(np.mean((y - x_clean)**2)):.3f} -> denoised {np.sqrt(np.mean((x - x_clean)**2)):.3f}")

The root-mean-square error drops from 0.507 to 0.123. ADMM converges faster per iteration but solves a linear system at each step. Chambolle–Pock only multiplies by D and D^\top: on a large two-dimensional image with general boundary conditions, it is the method whose iterations stay cheap and parallel.

5. Majorization–Minimization methods

Instead of minimizing a hard function directly, Majorization–Minimization (MM) minimizes, at each step, a simpler function that lies above it and touches it at the current point. Gradient descent, half-quadratic algorithms, iteratively reweighted least squares and EM are all instances.

The MM principle

Monotonicity holds without convexity, which makes MM popular for non-convex models. The art is to choose a majorant that is tight enough to make progress and simple enough to minimize.

Building majorants

The hyperbolic potential \varphi(t) = \sqrt{t^2 + \delta^2}, a smooth version of |t| that preserves edges, satisfies item 3 with \omega(s) = 1/\sqrt{s^2 + \delta^2}. So does the Cauchy potential \log(1 + t^2/\delta^2), which is not convex.

Item 3 is the majorant behind the multiplicative half-quadratic form of Geman and Reynolds (1992). The additive form of Geman and Yang (1995) corresponds instead to a majorant of constant curvature, as in item 1.

For a penalty \sum_i \varphi([Dx]_i), item 3 replaces each term by a weighted square, and the MM step solves a weighted least-squares problem: this is iteratively reweighted least squares (IRLS).

Convergence of quadratic MM

If moreover \nabla F is Lipschitz and F has the Kurdyka–Łojasiewicz property, which semi-algebraic functions satisfy, every bounded MM sequence converges as a whole, even without convexity (Attouch, Bolte and Svaiter, 2013); see ch. 9.

MM at scale: the memory-gradient subspace method

For a large image, the half-quadratic step solves a linear system with the matrix A(x_k) at every iteration, which is expensive. The subspace method minimizes the same majorant only over a small subspace.

The second direction acts as a memory, like the momentum of Algorithm 1.5, but its weight is optimized at each step rather than prescribed.

A separable majorant, with A(y) diagonal, is the other road to scale: Proposition 1.9 builds one from any symmetric matrix. Then every coordinate is updated independently, in parallel, which is how classical multiplicative algorithms such as the non-negative matrix factorization updates of Lee and Seung (2001) work.

import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import spsolve

rng = np.random.default_rng(0)
N, lam, delta = 2000, 1.0, 0.05
x_clean = np.repeat([0.0, 2.0, -1.0, 1.0, 0.0], N // 5)
y = x_clean + 0.5 * rng.standard_normal(N)
D = sp.diags([-np.ones(N - 1), np.ones(N - 1)], [0, 1], shape=(N - 1, N)).tocsr()

phi = lambda t: np.sqrt(t ** 2 + delta ** 2)          # hyperbolic (smoothed TV) potential
F = lambda x: 0.5 * np.sum((x - y) ** 2) + lam * phi(D @ x).sum()
grad = lambda x: x - y + lam * (D.T @ ((D @ x) / phi(D @ x)))
weights = lambda x: 1.0 / phi(D @ x)                  # omega(t) = phi'(t) / t

def mm_halfquadratic(iters):                           # exact minimisation of the majorant
    x = y.copy()
    for _ in range(iters):
        A = sp.eye(N) + lam * D.T @ sp.diags(weights(x)) @ D
        x = spsolve(A.tocsc(), y)
    return x

def mm_subspace(iters):                                # 3MG: majorant minimised on span(-g, x - x_prev)
    x, x_prev = y.copy(), y.copy()
    for k in range(iters):
        g, w = grad(x), weights(x)
        S = np.column_stack([-g, x - x_prev]) if k else -g[:, None]
        AS = S + lam * (D.T @ (w[:, None] * (D @ S)))  # A(x) S without forming A(x)
        u = np.linalg.solve(S.T @ AS, -S.T @ g)
        x_prev, x = x, x + S @ u
    return x

def gradient_descent(iters):
    L = 1 + 4 * lam / delta                            # Lipschitz bound: ||D||^2 <= 4, phi'' <= 1/delta
    x = y.copy()
    for _ in range(iters):
        x = x - grad(x) / L
    return x

F_star = F(mm_halfquadratic(300))
gap = lambda x: abs(F(x) - F_star)
for k in (10, 50, 200):
    print(f"k={k:4d}  MM-HQ {gap(mm_halfquadratic(k)):.2e}   "
          f"3MG {gap(mm_subspace(k)):.2e}   GD {gap(gradient_descent(k)):.2e}")
vals = []
x = y.copy()
for _ in range(10):                                    # MM is monotone: F decreases at every step
    A = sp.eye(N) + lam * D.T @ sp.diags(weights(x)) @ D
    x = spsolve(A.tocsc(), y)
    vals.append(F(x))
print("monotone decrease:", bool(np.all(np.diff(vals) <= 1e-12)))

The exact MM step converges in fewer iterations, but each one solves a sparse system. 3MG needs a few more iterations at the cost of a gradient step, and it is the version that scales to large two-dimensional problems. Gradient descent pays for the global constant L = 81, set by the sharpest edge, while the majorant adapts its curvature locally.

6. Stochastic and block-coordinate methods for large data sets

When the objective is an average over n = 10^9 samples, or has 10^9 variables, even one exact gradient is too expensive. This chapter replaces it by a cheap unbiased estimate (stochastic methods) or by a partial derivative (coordinate methods), and measures what is lost and how to recover it.

The objective is a finite sum, as in Example 1.1: F(x) = \frac{1}{n} \sum_{i=1}^n f_i(x) . One full gradient costs n sample gradients: a pass over the data.

Stochastic gradient

The rate O(1/t) is slow compared with the linear rate of gradient descent, but each step costs 1/n of a pass. The noise is the limit: with a constant step, SGD only reaches a neighborhood of x^\star whose size is proportional to the step and to the variance of the sample gradients.

The B sample gradients of a mini-batch are independent computations: on B workers, a mini-batch step takes the time of one sample gradient. This is the root of data parallelism (ch. 7).

Variance reduction

The variance of SGD does not vanish at the solution, because the individual gradients \nabla f_i(x^\star) are not zero. Variance-reduced methods correct each sample gradient with a reference point.

The direction is still unbiased, and its variance tends to zero as both x_t and \tilde x approach x^\star.

With \eta = 0.1/L and m = 50 L/\mu, the factor is \rho = \frac{1}{5 \times 0.8} + \frac{0.2}{0.8} = 0.5.

Each outer step costs n + 2m sample gradients. Reaching accuracy \varepsilon therefore costs O((n + L/\mu) \log(1/\varepsilon)) sample gradients, against O(n\, \kappa \log(1/\varepsilon)) for gradient descent: the condition number and the number of samples now add instead of multiplying.

SAGA (Defazio, Bach and Lacoste-Julien, 2014) keeps a table of the last gradient of each f_i instead of a snapshot. It needs no full passes and handles a proximal term.

import numpy as np
from scipy.optimize import minimize

rng = np.random.default_rng(0)
n, d, mu = 5000, 20, 1e-3
A = rng.standard_normal((n, d))
b = np.sign(A @ rng.standard_normal(d) + 0.5 * rng.standard_normal(n))

def F(x):
    return np.mean(np.logaddexp(0, -b * (A @ x))) + 0.5 * mu * x @ x
def grad_i(x, i):                                  # gradient of the i-th term f_i
    s = -b[i] / (1 + np.exp(b[i] * (A[i] @ x)))
    return s * A[i] + mu * x
def full_grad(x):
    s = -b / (1 + np.exp(b * (A @ x)))
    return A.T @ s / n + mu * x

F_star = minimize(F, np.zeros(d), jac=full_grad, method="L-BFGS-B", tol=1e-14).fun
L = np.linalg.norm(A, 2) ** 2 / (4 * n) + mu        # smoothness of F
L_max = np.max(np.sum(A ** 2, axis=1)) / 4 + mu     # smoothness of the worst f_i
passes = 20                                          # budget: 20 passes over the data

x = np.zeros(d)                                      # full gradient descent: 1 pass per step
for _ in range(passes):
    x -= full_grad(x) / L
gd = F(x) - F_star

x = np.zeros(d)                                      # SGD, decreasing step 1 / (L_max + mu t)
for t, i in enumerate(rng.integers(n, size=passes * n)):
    x -= grad_i(x, i) / (L_max + mu * t)
sgd = F(x) - F_star

eta = 0.5 / L_max                                    # SVRG: snapshot every n inner steps
x_tilde, used = np.zeros(d), 0
while used + 3 <= passes:                            # each outer loop costs 3n sample gradients
    g_tilde, x = full_grad(x_tilde), x_tilde.copy()
    for i in rng.integers(n, size=n):
        x -= eta * (grad_i(x, i) - grad_i(x_tilde, i) + g_tilde)
    x_tilde, used = x, used + 3
svrg = F(x_tilde) - F_star
print(f"after {passes} passes: GD {gd:.1e}   SGD {sgd:.1e}   SVRG {svrg:.1e}")
print(f"condition number L/mu = {L / mu:.0f},  L_max/L = {L_max / L:.1f}")

After 20 passes the objective gaps are \text{GD: } 2.5 \times 10^{-2}, \qquad \text{SGD: } 2.0 \times 10^{-3}, \qquad \text{SVRG: } 3.0 \times 10^{-8} . The condition number is 285 and L_{\max}/L = 41.4. SGD makes fast early progress and then stalls on its variance; SVRG keeps the cheap steps and converges linearly.

Randomized block-coordinate descent

When the number of variables is huge, one can update a single block of coordinates at a time. Write x = (x_1, \dots, x_p) in blocks and let L_i be the Lipschitz constant of the block gradient \nabla_i F with respect to changes of block i.

One block step costs about 1/p of a gradient when partial derivatives are cheap, for instance for F(x) = f(Ax) with the residual kept up to date. Per pass, the contraction is about e^{-\mu/L_{\max}}, against e^{-\mu/L} for gradient descent.

Since L_{\max} \leq L \leq p L_{\max}, the gain per pass lies between 1 and p. There is no gain for a separable F, where L = L_{\max}. The gain is largest when the coordinates are strongly coupled, for instance for least squares with highly correlated columns, where L approaches p L_{\max}.

Nesterov (2012) gives the general analysis with non-uniform sampling; the proximal and parallel versions (Richtárik and Takáč, 2014 and 2016) lead to ch. 7.

7. Parallel optimization

A multicore machine or a cluster offers p processors. This chapter asks which parts of an optimization algorithm can run simultaneously, how much time that saves, and what must change in the algorithm when updates happen at the same time.

Architectures and speedup

With s = 0.05 and p = 32, the speedup is only 1/(0.05 + 0.95/32) = 12.5, and no number of processors exceeds 20. The sequential part of an optimization algorithm is typically a reduction (an average, a norm, a step-size computation) or a synchronization.

Product space and the parallel proximal algorithm

Many problems are sums \sum_{i=1}^m f_i(x) of terms with simple proxes, while the prox of the sum is intractable. The idea is to give each term its own copy of the variable and to force the copies to agree.

Each iteration costs one prox per term, all simultaneous, and a single reduction of N numbers: one all-reduce in distributed memory (Proposition 7.3).

import numpy as np

rng = np.random.default_rng(0)
N, lam = 256, 0.05
x_true = np.zeros(N)
x_true[rng.choice(N, 12, replace=False)] = rng.uniform(0.5, 2.0, 12)   # positive spikes
h = np.exp(-0.5 * (np.arange(-8, 9) / 2.5) ** 2); h /= h.sum()         # Gaussian blur
Hf = np.fft.fft(np.roll(np.pad(h, (0, N - h.size)), -8))               # circular convolution
H = lambda v: np.real(np.fft.ifft(Hf * np.fft.fft(v)))
Ht = lambda v: np.real(np.fft.ifft(np.conj(Hf) * np.fft.fft(v)))
y = H(x_true) + 0.01 * rng.standard_normal(N)
F = lambda x: 0.5 * np.sum((H(x) - y) ** 2) + lam * np.abs(x).sum()   # on x >= 0

# the three proximity operators, each computable on its own core
def prox_data(v, g):        # (I + g H^T H)^{-1}(v + g H^T y), diagonal in Fourier
    return np.real(np.fft.ifft(np.fft.fft(v + g * Ht(y)) / (1 + g * np.abs(Hf) ** 2)))
prox_l1 = lambda v, g: np.sign(v) * np.maximum(np.abs(v) - g * lam, 0.0)
prox_pos = lambda v, g: np.maximum(v, 0.0)
proxes, w = [prox_data, prox_l1, prox_pos], np.full(3, 1 / 3)

def ppxa(gamma, iters=1000):            # Combettes and Pesquet (2008), lambda = 1
    z, x = [np.zeros(N) for _ in proxes], np.zeros(N)
    for _ in range(iters):
        p_i = [prox(zi, gamma / wi) for prox, zi, wi in zip(proxes, z, w)]  # parallel step
        p = sum(wi * pi for wi, pi in zip(w, p_i))                         # reduction
        z = [zi + (2 * p - x - pi) for zi, pi in zip(z, p_i)]
        x = p
    return x

v, t, u = np.zeros(N), 1.0, np.zeros(N)            # reference: FISTA with prox of l1 + positivity
for _ in range(20000):
    v_new = np.maximum(u - Ht(H(u) - y) - lam, 0.0)                   # here L = max |Hf|^2 = 1
    t_new = (1 + np.sqrt(1 + 4 * t * t)) / 2
    u, v, t = v_new + (t - 1) / t_new * (v_new - v), v_new, t_new
print(f"FISTA reference: F = {F(v):.6f}")
for gamma in (0.2, 1.0, 10.0):
    x = ppxa(gamma)
    print(f"PPXA gamma={gamma:5.1f}: F - F_ref = {F(x) - F(v):.1e}, "
          f"relative distance {np.linalg.norm(x - v) / np.linalg.norm(v):.1e}")

After 1000 iterations, the relative distance to the reference solution depends strongly on \gamma:

Theorem 7.6 holds for every \gamma > 0, but the speed does not: tuning \gamma to the scale of the problem is part of using any Douglas–Rachford-type method.

Parallel coordinate descent

Coordinate descent (Algorithm 6.7) is sequential: each step uses the latest iterate. Updating \tau blocks at once, one per core, is natural, but the simultaneous steps are all computed from the same point and can overshoot when the blocks interact.

The speedup \tau / \beta depends on the coupling:

• with p = 10^6 blocks, \tau = 32 cores and sparse coupling \omega = 10, \beta \approx 1.0003: the speedup is almost 32;
• with a dense problem, \omega = p, \beta = \tau: no speedup at all;
• with \omega = p/2, \beta \approx 16.5 and the speedup is about 2.

Exercise 12 shows on a small example why \beta = 1 can diverge.

Data parallelism and asynchrony

In data parallelism, each of p workers holds a shard of the data, computes the gradient of its part of the loss, and an all-reduce averages the gradients. With equal shards, this is exactly one mini-batch step (Proposition 6.3): the theory of ch. 6 applies unchanged, and only the communication of Proposition 7.3 is added.

import torch

torch.set_num_threads(1)
torch.manual_seed(0)
n, d, workers = 4000, 30, 4
X = torch.randn(n, d)
y = (X @ torch.randn(d) + 0.3 * torch.randn(n) > 0).float()
loss_fn = torch.nn.BCEWithLogitsLoss()

def train(n_workers, steps=200, lr=0.5):
    model = torch.nn.Linear(d, 1)
    torch.nn.init.zeros_(model.weight); torch.nn.init.zeros_(model.bias)
    shards = torch.arange(n).chunk(n_workers)         # data parallelism: one shard per worker
    for _ in range(steps):
        grads = []
        for idx in shards:                            # in practice: one process / GPU each
            model.zero_grad()
            loss_fn(model(X[idx]).squeeze(1), y[idx]).backward()
            grads.append([p.grad.clone() for p in model.parameters()])
        with torch.no_grad():                         # all-reduce: average the gradients
            for j, p in enumerate(model.parameters()):
                p -= lr * sum(g[j] for g in grads) / n_workers
    return torch.cat([p.detach().flatten() for p in model.parameters()])

w1, w4 = train(1), train(workers)
print(f"max |w(1 worker) - w({workers} workers)| = {(w1 - w4).abs().max():.1e}")

After 200 steps the two parameter vectors differ by 2.4 \times 10^{-7}, the rounding error of single precision. In a real run, the inner loop is replaced by one process per GPU and the average by an all-reduce, for instance with torch.distributed.

Asynchronous methods remove the barrier. Hogwild! (Recht et al., 2011) lets threads update a shared iterate without locks; it converges for sparse problems where conflicting writes are rare.

ARock (Peng et al., 2016) runs asynchronous block updates of a Krasnosel’skii–Mann iteration (Theorem 4.3) with delays of at most \tau iterations; here \tau is a delay, not the number of blocks of Algorithm 7.8. Its step size must shrink with the delay, roughly in proportion to 1/(1 + 2\tau/\sqrt{p}) for p blocks: asynchrony is almost free when \tau is small compared with \sqrt{p}. Exercise 15 computes the exact stability limit of a delayed gradient step on a quadratic.

8. Distributed optimization strategies

Now the data are split among m agents (machines, hospitals, phones) that cannot or must not pool them, because of size, privacy or regulation. Agent i holds a private function f_i, and together they must solve \min_{x \in \mathbb{R}^N} \; \sum_{i=1}^m f_i(x) . Communication, not computation, is the bottleneck. Two topologies are common: a star, where agents talk to a server, and a graph, where they only talk to their neighbours.

Star topology: consensus ADMM

Copy the variable on every agent and add the consensus constraint x_i = z. ADMM (Algorithm 4.6) on this formulation gives a method where each agent solves a local problem and the server only averages.

Convergence is Theorem 4.7, with f(x_1, \dots, x_m) = \sum_i f_i(x_i), g = 0 and the constraint x_i - z = 0. Each round costs one local prox per agent and one exchange of N numbers with the server. The data never leave the agents.

Graph topology: gossip

Decentralized gradient descent and its bias

The disagreement is penalized with weight 1/\alpha, so the bias shrinks with the step (Yuan, Ling and Yin, 2016). Decreasing steps remove it, but slowly.

Gradient tracking

With strongly convex smooth f_i and a small enough step, gradient tracking converges linearly to the exact solution (Nedić, Olshevsky and Shi, 2017). EXTRA (Shi et al., 2015) reaches the same goal with a correction based on two successive mixing matrices.

import numpy as np
import networkx as nx

rng = np.random.default_rng(0)
m, d, n_loc = 10, 5, 40                         # 10 agents, 40 private samples each
G = nx.cycle_graph(m)                           # communication graph: a ring
A = [rng.standard_normal((n_loc, d)) for i in range(m)]
b = [Ai @ rng.normal(1, 1, d) + 0.1 * rng.standard_normal(n_loc) for Ai in A]  # heterogeneous
grad = lambda i, x: A[i].T @ (A[i] @ x - b[i])  # f_i(x) = 0.5 ||A_i x - b_i||^2
x_star = np.linalg.lstsq(np.vstack(A), np.concatenate(b), rcond=None)[0]

W = np.zeros((m, m))                            # Metropolis-Hastings mixing matrix
for i, j in G.edges:
    W[i, j] = W[j, i] = 1 / (1 + max(G.degree[i], G.degree[j]))
W += np.diag(1 - W.sum(axis=1))
lam2 = np.sort(np.abs(np.linalg.eigvals(W)))[-2]
print(f"second largest |eigenvalue| of W: {lam2:.3f}")

L = max(np.linalg.norm(Ai, 2) ** 2 for Ai in A)
T = 3000
for c in (0.2, 0.1, 0.05):                      # DGD: mix, then local gradient step
    X, alpha = np.zeros((m, d)), c / L
    for _ in range(T):
        X = W @ X - alpha * np.array([grad(i, X[i]) for i in range(m)])
    print(f"DGD, alpha={c}/L:      max_i ||x_i - x*|| = {np.linalg.norm(X - x_star, axis=1).max():.2f}")

alpha = 0.1 / L
X = np.zeros((m, d))                            # gradient tracking
g_old = np.array([grad(i, X[i]) for i in range(m)]); S = g_old.copy()
for _ in range(T):
    X = W @ X - alpha * S
    g_new = np.array([grad(i, X[i]) for i in range(m)])
    S, g_old = W @ S + g_new - g_old, g_new
print(f"gradient tracking, 0.1/L: max_i ||x_i - x*|| = {np.linalg.norm(X - x_star, axis=1).max():.1e}")

rho = 50.0                                      # consensus ADMM with a central averaging step
X, U, z = np.zeros((m, d)), np.zeros((m, d)), np.zeros(d)
solves = [np.linalg.inv(Ai.T @ Ai + rho * np.eye(d)) for Ai in A]
for _ in range(50):
    X = np.array([solves[i] @ (A[i].T @ b[i] + rho * (z - U[i])) for i in range(m)])
    z = (X + U).mean(axis=0)
    U += X - z
print(f"consensus ADMM, 50 rounds:  ||z - x*|| = {np.linalg.norm(z - x_star):.1e}")

The mixing matrix has \lambda = 0.873, as in Example 8.4. The bias of DGD shrinks with the step, as Proposition 8.6 predicts, while gradient tracking is exact. Consensus ADMM is fastest, but it needs a server and exact local solves.

Federated learning

In federated averaging (McMahan et al., 2017), each client runs several local SGD steps before the server averages the models, to save communication rounds. With heterogeneous clients, the local steps drift toward the local minimizers, and the fixed point is biased: Exercise 14 computes the bias on two quadratics.

SCAFFOLD (Karimireddy et al., 2020) corrects the drift with control variates, much as gradient tracking corrects DGD. On the theory side, Scaman et al. (2017) characterize the optimal number of communication rounds for smooth strongly convex problems on a network: it grows with the square root of the condition number and with the inverse square root of the spectral gap of the gossip matrix.

9. Research directions and open problems

This self-study chapter surveys active questions around the course. It names only work whose authors and year are given in the references, and it does not try to be exhaustive.

Beyond convexity

Many imaging and learning models are non-convex: sparsity-inducing penalties such as \ell_0 or Cauchy potentials, matrix factorizations, neural networks. The Kurdyka–Łojasiewicz (KL) inequality has become the standard tool to prove that a whole sequence of iterates converges, not only that its gradients vanish.

• Attouch, Bolte and Svaiter (2013) give an abstract framework: sufficient decrease, a relative error condition and the KL property imply convergence of descent methods, including forward–backward.
• Bolte, Sabach and Teboulle (2014) build PALM, an alternating proximal linearized method for non-convex, non-smooth problems in two blocks of variables.
• Chouzenoux, Pesquet and Repetti (2014, 2016) extend forward–backward to variable metrics derived from MM majorants and to block-coordinate updates, under the KL property.

Open question: beyond convergence, which non-convex structures give explicit rates or guarantee convergence to a global minimizer?

Stochastic, incremental and randomized splitting

The deterministic splitting methods of ch. 4 touch all the data at each iteration. Several lines of work randomize them.

• Mairal (2015) proposes incremental MM: each step updates the majorant of one term of a finite sum, which brings the MM principle to large data sets.
• Combettes and Pesquet (2015) prove the convergence of block-coordinate fixed-point iterations in which a random subset of blocks is updated at each step, which covers randomized versions of forward–backward and Douglas–Rachford.
• Pesquet and Repetti (2015) and Bianchi, Hachem and Iutzeler (2016) derive randomized and asynchronous primal–dual algorithms for distributed optimization, in which each agent wakes up at random times and only talks to its neighbours.

The survey of Sun, Babu and Palomar (2017) covers MM in signal processing, communications and machine learning; Komodakis and Pesquet (2015) and Condat et al. (2023) survey primal–dual and splitting methods. Davis and Yin (2017) add a third operator to Douglas–Rachford, so that a smooth term enters through its gradient only.

Learning the algorithm: unrolling

Gregor and LeCun (2010) observed that a fixed number of ISTA iterations is a neural network whose weights can be learned: LISTA reaches a given accuracy with far fewer layers than ISTA needs iterations. Monga, Li and Eldar (2021) survey this algorithm unrolling. Bertocchi et al. (2020) unfold a proximal interior-point method for image restoration, learning its step sizes and regularization parameters.

Open questions: what guarantees (stability, robustness to the forward model, convergence as the depth grows) can an unrolled network keep from the algorithm it comes from?

Communication-efficient distributed learning

• Lian et al. (2017) show that decentralized SGD can match centralized mini-batch SGD while avoiding the bottleneck of a server.
• Koloskova et al. (2020) give a unified analysis of decentralized SGD with changing topologies and local updates.
• Alistarh et al. (2017) quantize gradients (QSGD) to cut the number of bits sent, with convergence guarantees.
• Mishchenko et al. (2022) show with ProxSkip that communicating only at random steps, with a control variate on the local updates, reduces the number of communication rounds from order \kappa to order \sqrt{\kappa} (times a logarithm) in the strongly convex case.
• Kairouz et al. (2021) collect the open problems of federated learning: heterogeneity of clients, privacy, robustness to faulty or adversarial clients, fairness.

Open problems

1. Tuning without knowledge of the problem. Douglas–Rachford, ADMM and PPXA converge for every step \gamma > 0, but their speed varies by orders of magnitude (Example 7.7). Adaptive and parameter-free rules with guarantees remain an active topic.
2. Heterogeneity versus communication. The bias of local steps (Exercise 14) and of decentralized gradients (Proposition 8.6) can be corrected, but the optimal trade-off between local computation and communication under strong heterogeneity and non-convexity is not settled.
3. Asynchrony with non-smooth and non-convex terms. Bounded-delay analyses exist for convex fixed-point iterations; extending them to the non-convex settings of deep learning, with realistic delays, is open.
4. MM for very large models. Majorants with cheap minimizers are well understood for imaging; building them for deep networks, where curvature varies widely between layers, is a natural but difficult direction.

Exercises

Levels: ★ apply the course, ★★ combine results or prove a short statement, ★★★ open-ended or research-flavoured. Exercises 6, 13, 14 and 15 are computational.

1. ★ Proximity operators by hand. Compute:

   • 1. the prox of |\cdot| (\gamma = 1) at v = 3, -0.4 and 1.2;
   • 2. the prox of \gamma f for f(x) = \frac{a}{2} x^2 on \mathbb{R}, with a > 0;
   • 3. the prox of \gamma \|\cdot\|_2 at v = (3, 4), for \gamma = 1 and \gamma = 6;
   • 4. the projection of (2, -0.5, 0.3) onto [0, 1]^3.

   Solution: (a) soft thresholding gives 2, 0 and 0.2. (b) Setting the derivative of \frac{a}{2}x^2 + \frac{1}{2\gamma}(x - v)^2 to zero gives \mathrm{prox}_{\gamma f}(v) = \frac{v}{1 + \gamma a} .

   3. Here \|v\| = 5. For \gamma = 1 the result is (1 - \frac15)(3, 4) = (2.4, 3.2); for \gamma = 6 \geq 5 it is 0 (Example 2.3, item 4).

   4. Clipping each coordinate gives (1, 0, 0.3).

2. ★ Step sizes and iteration counts. Let f(x) = \frac12 x^\top Q x - c^\top x with Q = \mathrm{diag}(1, 100). Give L, \mu and \kappa. How many iterations of gradient descent with step 1/L guarantee, by Theorem 0.7, \|x_k - x^\star\|^2 \leq 10^{-6}\, \|x_0 - x^\star\|^2 \, ? Compare with the order (1 - 1/\sqrt{\kappa})^k of the accelerated method.

   Solution: L = 100, \mu = 1, \kappa = 100. Theorem 0.7 needs 0.99^k \leq 10^{-6}: k \geq \frac{\ln 10^6}{-\ln 0.99} = \frac{13.82}{0.01005} \approx 1375 . For the accelerated rate, 0.9^k \leq 10^{-6} gives k \geq 13.82 / 0.1054 \approx 132: ten times fewer, the factor \sqrt{\kappa}.

3. ★ The Huber function. Derive the Moreau envelope of |\cdot| (Example 2.8) from Proposition 2.2, and check the gradient formula of Proposition 2.7. Evaluate the envelope and its gradient for \gamma = 0.5 at v = 0.2 and v = 2.

   Solution: for |v| \leq \gamma the prox is 0 and the value v^2/(2\gamma). Otherwise the prox is v - \gamma\,\mathrm{sign}(v) and the value |v| - \gamma/2.

   The gradient (v - \mathrm{prox}(v))/\gamma is v/\gamma in the first case and \mathrm{sign}(v) in the second: it is continuous and \frac1\gamma-Lipschitz.

   With \gamma = 0.5, at v = 0.2 the value is 0.04 and the gradient 0.4. At v = 2 the value is 1.75 and the gradient 1.

4. ★ One forward–backward step by hand. In the lasso, take A = \mathrm{diag}(2, 1), b = (4, 0.5) and \lambda = 1. Compute L, one ISTA step from x_0 = 0, and the objective before and after. Check optimality with Example 0.4.

   Solution: L = \|A\|^2 = 4. The gradient at 0 is -A^\top b = (-8, -0.5), so the forward step gives (2, 0.125). Soft thresholding at \lambda / L = 0.25 gives x_1 = (1.75,\ 0) . The objective goes from 8.125 to 0.5\,(0.25 + 0.25) + 1.75 = 2.

   Optimality: A^\top(b - Ax_1) = (1, 0.5). The first coordinate equals \lambda\,\mathrm{sign}(1.75), and the second satisfies |0.5| \leq \lambda where the solution is zero. So x_1 is optimal: one step suffices because the first coordinate has curvature exactly L.

5. ★ Gossip on a path of three agents. Write the Metropolis–Hastings matrix of the path 1 - 2 - 3, its eigenvalues and \lambda. Starting from x^0 = (3, 0, 0), compute two gossip rounds and the error norms. How many rounds reduce the error by 10^6?

   Solution: the degrees are 1, 2, 1, every edge weight is 1/3, and W = \begin{pmatrix} 2/3 & 1/3 & 0 \\ 1/3 & 1/3 & 1/3 \\ 0 & 1/3 & 2/3 \end{pmatrix} . The eigenvalues are 1, then 2/3 with eigenvector (1, 0, -1), then 0 with eigenvector (1, -2, 1), so \lambda = 2/3.

   The rounds give x^1 = (2, 1, 0) and x^2 = (5/3, 1, 1/3). The error norms are \sqrt6 \approx 2.45, \sqrt2 \approx 1.41 and \sqrt{8}/3 \approx 0.94.

   From round 1 on, the error is an eigenvector and shrinks by exactly 2/3 per round. A reduction by 10^6 needs 13.82/\ln 1.5 \approx 35 rounds.

6. ★ Randomized coordinate descent for the lasso, in Python. On the lasso of Example 3.8, implement coordinate descent with soft thresholding (Algorithm 6.7 with g = \lambda\|\cdot\|_1), keeping the residual r = Ax - b up to date. After 1, 10 and 50 epochs of d coordinate steps, compare with the same number of ISTA iterations.

   Hints and expected result: the step on coordinate j is a soft thresholding at \lambda/L_j of x_j - a_j^\top r / L_j, with L_j = \|a_j\|^2 \approx 1. It is followed by the update r \leftarrow r + a_j\,\Delta x_j, and one epoch costs about one ISTA iteration.

   Drawing the coordinates with the generator rng that made the data, the gaps are about 1.9 and 0.47 after 1 and 10 epochs, against 2.7 and 0.52 for ISTA. After 50 epochs, coordinate descent is near 2 \times 10^{-5} while ISTA is still at 0.29: the coordinate constants are about 1, while L \approx 6.5 (Theorem 6.8).

   Other seeds give between 1.1 and 1.8 after one epoch, and anywhere between 10^{-9} and 6 \times 10^{-3} after 50 epochs: average several runs before drawing conclusions.

7. ★★ The prox of the Euclidean norm, by Moreau’s identity. Prove the formula of Example 2.3, item 4, with Theorem 2.6. Deduce the prox of a group-lasso penalty \sum_g \|x_g\|_2 over disjoint groups.

   Solution: the conjugate of \|\cdot\|_2 is the indicator of the unit Euclidean ball B (Example 0.9), and \gamma^{-1}\iota_B = \iota_B. Theorem 2.6 gives \mathrm{prox}_{\gamma \|\cdot\|_2}(v) = v - \gamma P_B(v/\gamma) . If \|v\| \leq \gamma, then P_B(v/\gamma) = v/\gamma and the prox is 0. Otherwise P_B(v/\gamma) = v/\|v\|, and the prox is (1 - \gamma/\|v\|)\,v.

   For disjoint groups the penalty is separable, and the prox applies this shrinkage group by group (Example 2.3, item 5).

8. ★★ Fixed points do not depend on the step, stability does. Take f(x) = \frac{L}{2}x^2 on \mathbb{R} and g = 0.

   • 1. Show that forward–backward diverges for every x_0 \neq 0 when \gamma > 2/L, although 0 is a fixed point for every \gamma.
   • 2. Show that the forward–backward operator is nonexpansive exactly when \gamma \leq 2/L, and averaged when \gamma < 2/L.

   Solution: (a) x_{k+1} = (1 - \gamma L)\, x_k, and |1 - \gamma L| > 1 when \gamma > 2/L. The fixed-point equation x = (1 - \gamma L)x has the unique solution 0 for every \gamma > 0 (Proposition 3.2).

   2. The operator is multiplication by 1 - \gamma L, nonexpansive when |1 - \gamma L| \leq 1. For \gamma < 2/L it equals (1 - \alpha) I + \alpha (-I) with \alpha = \gamma L / 2 < 1, as Proposition 4.2 predicts.

9. ★★ Iteratively reweighted least squares as MM. Derive the MM iteration given by the half-quadratic majorant (Proposition 5.3, item 3) for F(x) = \tfrac12\|Ax - b\|^2 + \lambda \sum_i \sqrt{x_i^2 + \delta^2} . Run it by hand in the scalar case A = 1, b = 2, \lambda = 1, \delta = 0.1, from x_0 = 2, and compare the limit with the lasso solution (\delta = 0).

   Solution: with \omega(s) = 1/\sqrt{s^2 + \delta^2} and W_k = \mathrm{diag}(\omega(x_{k,i})), the majorant is a quadratic in x, minimized by x_{k+1} = \big(A^\top A + \lambda W_k\big)^{-1} A^\top b . In the scalar case, x_{k+1} = 2 / (1 + 1/\sqrt{x_k^2 + 0.01}). The iterates are 1.334, 1.144, 1.069, 1.036, and they decrease monotonically to 1.0049. The lasso solution is the soft thresholding of 2 at 1, that is 1.

   As \delta \to 0, the weights 1/|x_i| blow up at zero coordinates: an exact zero can never leave zero, which is why \delta > 0 is kept.

10. ★★ ADMM for the lasso at scale. Derive the three updates of Example 4.8 from Algorithm 4.6. When A is n \times d with n \ll d, show how to perform the x-update with one factorization of an n \times n matrix, and give the cost per iteration.

    Solution: Algorithm 4.6 with f(x) = \frac12\|Ax - b\|^2, g = \lambda\|\cdot\|_1 and the constraint x - z = 0 gives the x-step of Example 4.8.

    The z-step is soft thresholding of x_{k+1} + u_k at \lambda/\rho, and u accumulates x_{k+1} - z_{k+1}.

    The matrix inversion lemma gives (A^\top A + \rho I)^{-1} = \frac1\rho \Big( I - A^\top (\rho I + A A^\top)^{-1} A \Big) . Factorizing \rho I + AA^\top costs O(n^2 d) once. Each iteration then costs O(nd): two products by A or A^\top, and two triangular solves of size n.

11. ★★ The bias of DGD on two agents. Let f_1(x) = \frac12 (x - a)^2 and f_2(x) = \frac12 (x + a)^2, so that x^\star = 0, and let W have off-diagonal entries w \in (0, \frac12].

    Find the fixed point of DGD with step \alpha and check it against Proposition 8.6. Evaluate it for a = 1, w = 0.25 and \alpha = 0.1.

    Solution: by symmetry, the fixed point is (s, -s). The first equation reads s = (1 - w)s - ws - \alpha(s - a), hence s = \frac{\alpha a}{2w + \alpha} . Minimizing the penalized function of Proposition 8.6 gives the same point: J_\alpha(x_1, x_2) = \tfrac12(x_1 - a)^2 + \tfrac12(x_2 + a)^2 + \tfrac{w}{2\alpha}(x_1 - x_2)^2 . With the numbers, s = 0.1/0.6 = 1/6: each agent is off by 0.17. For small \alpha, the error is of order \alpha a / (2w).

12. ★★ Why parallel coordinate steps need damping. Let F(x) = \frac12 (\mathbf{1}^\top x)^2 on \mathbb{R}^p, so that every coordinate constant is L_i = 1, and write s = \mathbf{1}^\top x.

    • 1. Run Algorithm 7.8 with \tau coordinates and \beta = 1, and show that s^+ = (1 - \tau)s.
    • 2. Compute \beta from Theorem 7.9 and show that the method then converges in one step.
    • 3. For p = 1000, \tau = 32 and \omega = 4, compute \beta and the speedup.

    Solution: (a) \nabla_i F(x) = s for every i, so each updated coordinate loses s/\beta, and s^+ = (1 - \tau/\beta)s.

    With \beta = 1 and \tau \geq 3, the factor satisfies |1 - \tau| > 1: the method diverges.

    2. Every term couples all coordinates, so \omega = p and \beta = 1 + (p - 1)(\tau - 1)/(p - 1) = \tau. Then s^+ = 0: a minimizer.

    3. \beta = 1 + 3 \times 31/999 \approx 1.093, and the speedup is 32/1.093 \approx 29.

13. ★★ Spectral gap of rings, in Python. With NetworkX, build the Metropolis–Hastings matrix of rings with m = 10, 20, 40, 80 agents, compute 1 - \lambda and compare with 4\pi^2/(3m^2). Repeat for a \sqrt m \times \sqrt m torus and a complete graph, and estimate how the number of gossip rounds of Example 8.4 scales with m.

    Hints and expected result: for the rings, 1 - \lambda is 0.1273, 0.0326, 0.0082 and 0.0021, against the approximations 0.1316, 0.0329, 0.0082 and 0.0021. The rounds grow like m^2 on a ring.

    For a torus, 1 - \lambda decays like 1/m, so the rounds grow like m.

    On a complete graph, the Metropolis weights are all 1/m, W is the averaging matrix and \lambda = 0. One round suffices, at the price of m - 1 messages per agent.

14. ★★★ Client drift in federated averaging. Two clients hold f_i(x) = \frac{c_i}{2}(x - a_i)^2. Each round, both start from the server value x and run K gradient steps of size \eta on their own f_i; the server then averages the results.

    • 1. Find the fixed point and compare it with x^\star = \sum_i c_i a_i / \sum_i c_i.
    • 2. Evaluate for c = (1, 4), a = (0, 1), \eta = 0.1, K = 10, and check by simulation.
    • 3. Propose and test a correction.

    Hints and expected result: (a) with q_i = (1 - \eta c_i)^K, client i returns a_i + q_i(x - a_i), and the fixed point is \hat x = \frac{\sum_i (1 - q_i)\, a_i}{\sum_i (1 - q_i)} . It equals x^\star for K = 1, where 1 - q_i = \eta c_i, but not in general.

    2. q = (0.349, 0.006), so \hat x = 0.604 while x^\star = 0.8; a simulation of 500 rounds converges to 0.6041.

    3. A control variate in the spirit of SCAFFOLD (Karimireddy et al., 2020) adds \nabla F(x) - \nabla f_i(x), computed at the start of the round, to every local gradient. The fixed point becomes exact; with the numbers of (b), the error is multiplied by about -0.12 per round, against 0.75 for plain gradient steps with K = 1.

15. ★★★ How much delay can a gradient step tolerate? Asynchronous workers apply gradients computed on stale iterates. On f(x) = \frac{L}{2}x^2, consider the delayed iteration x_{k+1} = x_k - \gamma L\, x_{k - \tau} .

    • 1. For \tau = 1, show that the iteration is stable exactly when \gamma < 1/L.
    • 2. Check numerically, with the roots of the characteristic polynomial, that the limit for general \tau is \gamma L < 2\sin\big(\pi / (2(2\tau + 1))\big), and relate it to the step-size rule of ARock (ch. 7).

    Hints and expected result: (a) the characteristic polynomial is z^2 - z + a with a = \gamma L. For 0 < a \leq \frac14 both roots are real and lie in (0, 1).

    For a > \frac14 they are complex conjugate with |z|^2 = a, so the iteration is stable if and only if a < 1: the delay halves the admissible step.

    2. Build the polynomial z^{\tau + 1} - z^{\tau} + a and compute the largest root modulus with numpy.roots, for a just below and above the bound. For \tau = 2 and \tau = 5, the bounds 0.618 and 0.285 are where it crosses 1.

    For large \tau the bound behaves like \pi/(2\tau): the step must shrink like 1/\tau, the behaviour that the ARock step size encodes.

Projects and assessment

The official assessment, per master-mva.com, combines two lab reports and an exam; no weights are displayed, and the Université Paris-Saclay course detail classifies it as “Evaluation Continue non Intégrale” (partial continuous assessment). The solo subjects below are sized like an extended lab report; the group subjects are longer studies for students who want to go further. The written exam proposal follows the 2 h format announced on master-mva.com. None of these proposals is official material.

Solo projects

1. [Solo · ≈ 20 h] FISTA for image deblurring: a paper study. Reproduce the main experiment of Beck and Teboulle (2009): ISTA against FISTA on a deblurring problem with \ell_1 regularization in an orthonormal transform.
   • Implement: a Gaussian blur by FFT, an orthonormal 2-D DCT (scipy.fft.dctn with norm="ortho") as sparsifying transform, ISTA, FISTA and FISTA with backtracking.
   • Data: a grey-level version of an image bundled with scikit-learn (sklearn.datasets.load_sample_image).
   • Success criterion: log–log plots of F(x_k) - F^\star with slopes close to -1 and -2 (Theorems 3.4 and 3.6), and the PSNR against the iteration number.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 2, 3.
2. [Solo · ≈ 20 h] Subspace Majorization–Minimization for image restoration: a paper study. Reimplement the 3MG algorithm of Chouzenoux, Idier and Moussaoui (2011) on two-dimensional deblurring with a hyperbolic edge-preserving penalty.
   • Compare: half-quadratic MM with conjugate-gradient inner solves, 3MG, nonlinear conjugate gradient and L-BFGS (scipy.optimize), all from the same start.
   • Success criterion: a table of time and iterations to reach a relative objective accuracy of 10^{-6}, and a check of monotone decrease (Proposition 5.2).
   • Deliverables: code, report, talk. Chapters 1, 5.
3. [Solo · ≈ 15 h] Variance reduction in practice. Compare SGD, SVRG (Johnson and Zhang, 2013) and SAGA (Defazio, Bach and Lacoste-Julien, 2014) on \ell_2-regularized logistic regression with n = 10^5 synthetic samples (sklearn.datasets.make_classification), and the proximal version of SAGA with an \ell_1 term.
   • Success criterion: suboptimality against the number of passes over 30 passes, a linear rate for SVRG and SAGA, and an experimental check of the condition \rho < 1 of Theorem 6.5.
   • Deliverables: code, report, talk. Chapters 3, 6.

Group projects

1. [Group of 2–3 · ≈ 30 h each] Distributed ADMM: a paper study. Reproduce the consensus and sharing formulations of Boyd et al. (2011) for \ell_1-regularized logistic regression, and compare them with decentralized methods.
   • Implementation: several processes on one machine with torch.distributed (gloo backend, CPU) or multiprocessing; synthetic data split across 2 to 8 workers, homogeneous or heterogeneous.
   • Work split: one member writes consensus ADMM with a server (Algorithm 8.1); one writes DGD and gradient tracking on a ring (Algorithms 8.5 and 8.7); one builds the communication layer, the timings and the comparison with Amdahl’s law (Proposition 7.2).
   • Success criterion: curves of accuracy against communication rounds and against wall-clock time, and the measured speedup for 1, 2, 4 and 8 workers, explained by the sequential fraction.
   • Deliverables: code, a 6-page report, a 15-minute talk. Chapters 4, 7, 8.
2. [Group of 2–3 · ≈ 30 h each] Parallel proximal algorithms for image restoration. Solve deblurring with three terms (data fidelity, total variation, a box constraint) with PPXA (Combettes and Pesquet, 2008), Chambolle–Pock or Condat–Vũ, and ADMM.
   • Data: images bundled with scikit-learn, blurred and corrupted with Gaussian noise.
   • Work split: one algorithm per member, with the proxes of PPXA computed in parallel processes (concurrent.futures), plus a shared harness for the objective, the PSNR and the timings.
   • Success criterion: all methods reach the same objective value to 10^{-6} relative accuracy; a study of the sensitivity to \gamma, \tau, \sigma and \rho; the time per iteration with 1 and 3 workers.
   • Deliverables: code, report, talk. Chapters 2, 4, 7.
3. [Group of 2–3 · ≈ 30 h each] Federated learning with heterogeneous clients: a paper study. Reproduce the qualitative findings of Karimireddy et al. (2020) comparing FedAvg (McMahan et al., 2017) and SCAFFOLD.
   • Data: the 8×8 digits of scikit-learn split across 10 clients by label, so that each client sees only two classes, with multinomial logistic regression or a small PyTorch network.
   • Work split: the simulation framework and FedAvg; SCAFFOLD; the experiments with the number of local steps K and the degree of heterogeneity, with the quadratic analysis of Exercise 14 as a sanity check.
   • Success criterion: curves of test accuracy against communication rounds that show the drift of FedAvg as K grows and its correction by SCAFFOLD.
   • Deliverables: code, report, talk. Chapters 6, 8.

Written exam (proposal)

Format: 2 hours, as announced on master-mva.com; one handwritten A4 sheet allowed, no electronic device; four independent questions.

1. Course question (5 points). Define the proximity operator of f \in \Gamma_0(\mathbb{R}^N) and prove that it is firmly nonexpansive. State the Moreau decomposition, and use it to compute the prox of \gamma \|\cdot\|_\infty in terms of a projection.

   Expected: Definition 2.1 with Proposition 2.2; the monotonicity argument of Proposition 2.5; Theorem 2.6. The conjugate of \|\cdot\|_\infty is the indicator of the unit \ell_1 ball B_1, so the prox of \gamma\|\cdot\|_\infty at v is v - \gamma P_{B_1}(v/\gamma).

2. Computation (6 points). Consider the problem \min_{x \in \mathbb{R}^2} \; \tfrac12 \|x - y\|^2 + \lambda |x_1 - x_2|, \qquad y = (3, 0), \quad \lambda = 1 .

   1. Solve it with the optimality conditions. (b) Write its Fenchel dual (Proposition 0.10) and solve it. (c) Write the Chambolle–Pock iteration and the condition on \tau, \sigma.

   Expected:

   • 1. With D = (1, -1), the guess x_1 > x_2 gives x - y + \lambda D^\top = 0, so x^\star = (2, 1), which is consistent. The optimal value is 1 + 1 = 2.
   • 2. The dual function is -u^2 + 3u on |u| \leq 1, maximized at u^\star = 1 with value 2: there is no duality gap, and x^\star = y - D^\top u^\star.
   • 3. Algorithm 4.9, with a clipping to [-1, 1] in the dual step and the prox of \frac12\|\cdot - y\|^2 in the primal one. Since \|D\|^2 = 2, the condition is \tau\sigma < 1/2.

3. Short proof (5 points). (a) Prove that decentralized gradient descent with a constant step is gradient descent on the penalized function J_\alpha of Proposition 8.6. (b) Deduce that its fixed point is not consensual on two agents with f_1(x) = \frac12(x - 1)^2 and f_2(x) = \frac12(x + 1)^2, and compute it for w = \alpha = 0.25.

   Expected: (a) the gradient computation of Proposition 8.6. (b) A consensual fixed point would need \nabla f_1(x) = \nabla f_2(x) = 0, impossible; Exercise 11 gives s = \alpha/(2w + \alpha) = 1/3, so the fixed point is (1/3, -1/3).

4. Problem (4 points). An iteration of PPXA with m = 8 terms costs 8 proxes of 10 ms each and an average of 2 ms that is sequential. (a) Give the time per iteration on 1 and 8 cores, and the speedup. (b) Check it with Amdahl’s law. (c) On a cluster, each iteration also needs a ring all-reduce of N = 10^7 doubles over 8 workers, on links of 10 GB/s; give the new time and speedup.

   Expected: (a) 82 ms and 12 ms, a speedup of 6.8. (b) s = 2/82 \approx 0.024 and 1/(s + (1 - s)/8) \approx 6.8. (c) By Proposition 7.3 each worker sends \frac{2 \times 7}{8} \times 10^7 doubles, that is 1.4 \times 10^8 bytes, about 14 ms; the iteration takes about 26 ms and the speedup falls to about 3.2.

Grading rubric (proposal)

References

Official references (master-mva.com)

• H. H. Bauschke and P. L. Combettes. Convex Analysis and Monotone Operator Theory in Hilbert Spaces. Springer, New York, 2011 (second edition 2017).
• N. Parikh and S. Boyd. Proximal algorithms. Foundations and Trends in Optimization, 1(3):127–239, 2014.

Textbooks and surveys

• A. Beck. First-Order Methods in Optimization. MOS-SIAM Series on Optimization, SIAM, 2017.
• D. P. Bertsekas and J. N. Tsitsiklis. Parallel and Distributed Computation: Numerical Methods. Prentice Hall, 1989.
• L. Bottou, F. E. Curtis and J. Nocedal. Optimization methods for large-scale machine learning. SIAM Review, 60(2):223–311, 2018.
• S. Boyd, N. Parikh, E. Chu, B. Peleato and J. Eckstein. Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends in Machine Learning, 3(1):1–122, 2011.
• S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press, 2004.
• P. L. Combettes and J.-C. Pesquet. Proximal splitting methods in signal processing. In H. H. Bauschke et al. (eds.), Fixed-Point Algorithms for Inverse Problems in Science and Engineering, Springer, 185–212, 2011.
• L. Condat, D. Kitahara, A. Contreras and A. Hirabayashi. Proximal splitting algorithms for convex optimization: a tour of recent advances, with new twists. SIAM Review, 65(2):375–435, 2023.
• P. Kairouz, H. B. McMahan et al. Advances and open problems in federated learning. Foundations and Trends in Machine Learning, 14(1–2):1–210, 2021.
• N. Komodakis and J.-C. Pesquet. Playing with duality: an overview of recent primal–dual approaches for solving large-scale optimization problems. IEEE Signal Processing Magazine, 32(6):31–54, 2015.
• V. Monga, Y. Li and Y. C. Eldar. Algorithm unrolling: interpretable, efficient deep learning for signal and image processing. IEEE Signal Processing Magazine, 38(2):18–44, 2021.
• Y. Nesterov. Introductory Lectures on Convex Optimization: A Basic Course. Kluwer Academic Publishers, 2004.
• Y. Sun, P. Babu and D. P. Palomar. Majorization-minimization algorithms in signal processing, communications, and machine learning. IEEE Transactions on Signal Processing, 65(3):794–816, 2017.

Proximal and splitting algorithms (ch. 2–4)

• A. Beck and M. Teboulle. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences, 2(1):183–202, 2009.
• A. Chambolle and C. Dossal. On the convergence of the iterates of the “fast iterative shrinkage/thresholding algorithm”. Journal of Optimization Theory and Applications, 166(3):968–982, 2015.
• A. Chambolle and T. Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision, 40(1):120–145, 2011.
• E. Chouzenoux, J.-C. Pesquet and A. Repetti. Variable metric forward–backward algorithm for minimizing the sum of a differentiable function and a convex function. Journal of Optimization Theory and Applications, 162(1):107–132, 2014.
• P. L. Combettes and V. R. Wajs. Signal recovery by proximal forward-backward splitting. Multiscale Modeling and Simulation, 4(4):1168–1200, 2005.
• L. Condat. A primal–dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms. Journal of Optimization Theory and Applications, 158(2):460–479, 2013.
• I. Daubechies, M. Defrise and C. De Mol. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Communications on Pure and Applied Mathematics, 57(11):1413–1457, 2004.
• D. Davis and W. Yin. A three-operator splitting scheme and its optimization applications. Set-Valued and Variational Analysis, 25(4):829–858, 2017.
• J. Eckstein and D. P. Bertsekas. On the Douglas–Rachford splitting method and the proximal point algorithm for maximal monotone operators. Mathematical Programming, 55:293–318, 1992.
• B. He and X. Yuan. Convergence analysis of primal-dual algorithms for a saddle-point problem: from contraction perspective. SIAM Journal on Imaging Sciences, 5(1):119–149, 2012.
• P.-L. Lions and B. Mercier. Splitting algorithms for the sum of two nonlinear operators. SIAM Journal on Numerical Analysis, 16(6):964–979, 1979.
• Y. Nesterov. A method for solving the convex programming problem with convergence rate O(1/k²). Soviet Mathematics Doklady, 27:372–376, 1983.
• R. T. Rockafellar. Monotone operators and the proximal point algorithm. SIAM Journal on Control and Optimization, 14(5):877–898, 1976.
• B. C. Vũ. A splitting algorithm for dual monotone inclusions involving cocoercive operators. Advances in Computational Mathematics, 38(3):667–681, 2013.

Majorization–Minimization (ch. 5)

• H. Attouch, J. Bolte and B. F. Svaiter. Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward–backward splitting, and regularized Gauss–Seidel methods. Mathematical Programming, 137(1–2):91–129, 2013.
• E. Chouzenoux, J. Idier and S. Moussaoui. A majorize–minimize strategy for subspace optimization applied to image restoration. IEEE Transactions on Image Processing, 20(6):1517–1528, 2011.
• D. Geman and G. Reynolds. Constrained restoration and the recovery of discontinuities. IEEE Transactions on Pattern Analysis and Machine Intelligence, 14(3):367–383, 1992.
• D. Geman and C. Yang. Nonlinear image recovery with half-quadratic regularization. IEEE Transactions on Image Processing, 4(7):932–946, 1995.
• D. D. Lee and H. S. Seung. Algorithms for non-negative matrix factorization. Advances in Neural Information Processing Systems 13, 2001.
• J. Mairal. Incremental majorization-minimization optimization with application to large-scale machine learning. SIAM Journal on Optimization, 25(2):829–855, 2015.

Stochastic, coordinate and parallel methods (ch. 6–7)

• P. L. Combettes and J.-C. Pesquet. A proximal decomposition method for solving convex variational inverse problems. Inverse Problems, 24(6):065014, 2008.
• A. Defazio, F. Bach and S. Lacoste-Julien. SAGA: a fast incremental gradient method with support for non-strongly convex composite objectives. Advances in Neural Information Processing Systems 27, 2014.
• R. Johnson and T. Zhang. Accelerating stochastic gradient descent using predictive variance reduction. Advances in Neural Information Processing Systems 26, 2013.
• Y. Nesterov. Efficiency of coordinate descent methods on huge-scale optimization problems. SIAM Journal on Optimization, 22(2):341–362, 2012.
• P. Patarasuk and X. Yuan. Bandwidth optimal all-reduce algorithms for clusters of workstations. Journal of Parallel and Distributed Computing, 69(2):117–124, 2009.
• Z. Peng, Y. Xu, M. Yan and W. Yin. ARock: an algorithmic framework for asynchronous parallel coordinate updates. SIAM Journal on Scientific Computing, 38(5):A2851–A2879, 2016.
• B. Recht, C. Ré, S. Wright and F. Niu. Hogwild!: a lock-free approach to parallelizing stochastic gradient descent. Advances in Neural Information Processing Systems 24, 2011.
• P. Richtárik and M. Takáč. Iteration complexity of randomized block-coordinate descent methods for minimizing a composite function. Mathematical Programming, 144(1–2):1–38, 2014.
• P. Richtárik and M. Takáč. Parallel coordinate descent methods for big data optimization. Mathematical Programming, 156(1–2):433–484, 2016.

Distributed and federated optimization (ch. 8)

• S. P. Karimireddy, S. Kale, M. Mohri, S. Reddi, S. Stich and A. T. Suresh. SCAFFOLD: stochastic controlled averaging for federated learning. Proceedings of ICML, 2020.
• B. McMahan, E. Moore, D. Ramage, S. Hampson and B. Agüera y Arcas. Communication-efficient learning of deep networks from decentralized data. Proceedings of AISTATS, 2017.
• A. Nedić and A. Ozdaglar. Distributed subgradient methods for multi-agent optimization. IEEE Transactions on Automatic Control, 54(1):48–61, 2009.
• A. Nedić, A. Olshevsky and W. Shi. Achieving geometric convergence for distributed optimization over time-varying graphs. SIAM Journal on Optimization, 27(4):2597–2633, 2017.
• G. Qu and N. Li. Harnessing smoothness to accelerate distributed optimization. IEEE Transactions on Control of Network Systems, 5(3):1245–1260, 2018.
• K. Scaman, F. Bach, S. Bubeck, Y. T. Lee and L. Massoulié. Optimal algorithms for smooth and strongly convex distributed optimization in networks. Proceedings of ICML, 2017.
• W. Shi, Q. Ling, G. Wu and W. Yin. EXTRA: an exact first-order algorithm for decentralized consensus optimization. SIAM Journal on Optimization, 25(2):944–966, 2015.
• K. Yuan, Q. Ling and W. Yin. On the convergence of decentralized gradient descent. SIAM Journal on Optimization, 26(3):1835–1854, 2016.

Research directions (ch. 9)

• D. Alistarh, D. Grubic, J. Li, R. Tomioka and M. Vojnovic. QSGD: communication-efficient SGD via gradient quantization and encoding. Advances in Neural Information Processing Systems 30, 2017.
• C. Bertocchi, E. Chouzenoux, M.-C. Corbineau, J.-C. Pesquet and M. Prato. Deep unfolding of a proximal interior point method for image restoration. Inverse Problems, 36(3):034005, 2020.
• P. Bianchi, W. Hachem and F. Iutzeler. A coordinate descent primal-dual algorithm and application to distributed asynchronous optimization. IEEE Transactions on Automatic Control, 61(10):2947–2957, 2016.
• J. Bolte, S. Sabach and M. Teboulle. Proximal alternating linearized minimization for nonconvex and nonsmooth problems. Mathematical Programming, 146(1–2):459–494, 2014.
• E. Chouzenoux, J.-C. Pesquet and A. Repetti. A block coordinate variable metric forward–backward algorithm. Journal of Global Optimization, 66(3):457–485, 2016.
• P. L. Combettes and J.-C. Pesquet. Stochastic quasi-Fejér block-coordinate fixed point iterations with random sweeping. SIAM Journal on Optimization, 25(2):1221–1248, 2015.
• K. Gregor and Y. LeCun. Learning fast approximations of sparse coding. Proceedings of ICML, 2010.
• A. Koloskova, N. Loizou, S. Boreiri, M. Jaggi and S. Stich. A unified theory of decentralized SGD with changing topology and local updates. Proceedings of ICML, 2020.
• X. Lian, C. Zhang, H. Zhang, C.-J. Hsieh, W. Zhang and J. Liu. Can decentralized algorithms outperform centralized algorithms? A case study for decentralized parallel stochastic gradient descent. Advances in Neural Information Processing Systems 30, 2017.
• K. Mishchenko, G. Malinovsky, S. Stich and P. Richtárik. ProxSkip: Yes! Local gradient steps provably lead to communication acceleration! Finally! Proceedings of ICML, 2022.
• J.-C. Pesquet and A. Repetti. A class of randomized primal-dual algorithms for distributed optimization. Journal of Nonlinear and Convex Analysis, 16(12):2453–2490, 2015.