QuantumCore

Lesson 32 · 100 min · core

Supervised Learning & Gradient Descent

Loading progress…

Supervised learning chooses a predictor by minimising the average loss on training data. This lesson formulates empirical-risk minimisation (ERM), proves that the training risk of an ERM solution is optimistically biased and decomposes the expected squared error into noise, bias2^2 and variance. It derives the gradients of the mean-squared-error loss of a linear model and of the binary cross-entropy of logistic regression coordinate by coordinate and in matrix form, proves the descent lemma for LL-smooth functions and the rate f(xk)−f∗≤L∥x0−x∗∥2/(2k)f(x_k)-f^\ast\le L\lVert x_0-x^\ast\rVert^2/(2k) of gradient descent with step 1/L1/L on convex objectives, shows that the iteration diverges on a quadratic once the step exceeds 2/λmax⁡2/\lambda_{\max}, and proves that minibatch gradients are unbiased with variance proportional to 1/B1/B — the same structure as the shot-noise-limited gradients of variational quantum circuits.

Objectives and key results

  • Formulate empirical-risk minimisation f^∈arg min⁡f∈F1m∑iℓ(f(xi),yi)\hat f\in\operatorname*{arg\,min}_{f\in\mathcal F}\tfrac1m\sum_i\ell(f(x_i),y_i), define the population risk and the generalisation gap, and prove that E[R^S(f^S)]≤inf⁡f∈FR(f)≤E[R(f^S)]\mathbb E[\hat R_S(\hat f_S)]\le\inf_{f\in\mathcal F}R(f)\le\mathbb E[R(\hat f_S)].

  • Derive the bias–variance decomposition E[(y−f^S(x))2]=σ2+bias2+variance\mathbb E[(y-\hat f_S(x))^2]=\sigma^2+\mathrm{bias}^2+\mathrm{variance} with every algebraic step.

  • Derive ∇wL=2m∑i(y^i−yi)xi=2mX⊤r\nabla_w\mathcal L=\tfrac2m\sum_i(\hat y_i-y_i)x_i=\tfrac2mX^\top r and ∂bL=2m∑i(y^i−yi)\partial_b\mathcal L=\tfrac2m\sum_i(\hat y_i-y_i) for the MSE loss, and ∇wL=1m∑i(σ(zi)−yi)xi\nabla_w\mathcal L=\tfrac1m\sum_i(\sigma(z_i)-y_i)x_i for the binary cross-entropy, including σ′=σ(1−σ)\sigma'=\sigma(1-\sigma).

  • Prove the descent lemma f(y)≤f(x)+∇f(x)⊤(y−x)+L2∥y−x∥2f(y)\le f(x)+\nabla f(x)^\top(y-x)+\tfrac L2\lVert y-x\rVert^2 and the O(1/k)O(1/k) convergence rate of gradient descent on convex LL-smooth functions with a full telescoping argument.

  • Show with the spectral theorem that gradient descent on a positive-definite quadratic converges for every start if and only if 0<η<2/λmax⁡0<\eta<2/\lambda_{\max}.

  • Prove that a minibatch gradient is an unbiased estimator of the full gradient with covariance Σ/B\Sigma/B, and compute the variance (1−μ2)/N(1-\mu^2)/N of an NN-shot Pauli expectation estimate.

  • Implement a vectorised MSE gradient in TypeScript and check it against finite differences.

definition · Empirical-risk minimisation
f^S∈arg min⁡f∈FR^S(f),R^S(f)=1m∑i=1mℓ(f(xi),yi),R(f)=E(x,y)∼D ℓ(f(x),y)\hat f_S\in\operatorname*{arg\,min}_{f\in\mathcal F}\hat R_S(f),\qquad \hat R_S(f)=\frac1m\sum_{i=1}^m\ell\big(f(x_i),y_i\big),\qquad R(f)=\mathbb E_{(x,y)\sim\mathcal D}\,\ell\big(f(x),y\big)
formula · Gradient of the mean-squared error
∇wL=2m∑i=1m(y^i−yi) xi=2mX⊤r,∂L∂b=2m∑i=1m(y^i−yi)\nabla_w\mathcal L=\frac2m\sum_{i=1}^m(\hat y_i-y_i)\,x_i=\frac2mX^\top r,\qquad \frac{\partial\mathcal L}{\partial b}=\frac2m\sum_{i=1}^m(\hat y_i-y_i)
formula · Gradient of the logistic (binary cross-entropy) loss
σ′(z)=σ(z)(1−σ(z)),∇wL=1m∑i=1m(σ(zi)−yi) xi\sigma'(z)=\sigma(z)\big(1-\sigma(z)\big),\qquad \nabla_w\mathcal L=\frac1m\sum_{i=1}^m\big(\sigma(z_i)-y_i\big)\,x_i
theorem · Descent lemma
f(y)≤f(x)+∇f(x)⊤(y−x)+L2∥y−x∥2,f(x−1L∇f(x))≤f(x)−12L∥∇f(x)∥2f(y)\le f(x)+\nabla f(x)^\top(y-x)+\frac L2\lVert y-x\rVert^2,\qquad f\Big(x-\tfrac1L\nabla f(x)\Big)\le f(x)-\frac1{2L}\lVert\nabla f(x)\rVert^2
theorem · O(1/k) rate of gradient descent
f convex, L-smooth, xk+1=xk−1L∇f(xk) ⟹ f(xk)−f∗≤L∥x0−x∗∥22kf\ \text{convex},\ L\text{-smooth},\ x_{k+1}=x_k-\tfrac1L\nabla f(x_k)\ \Longrightarrow\ f(x_k)-f^\ast\le\frac{L\lVert x_0-x^\ast\rVert^2}{2k}

Supervised learning and empirical-risk minimisation#

Quantum machine learning reuses the classical training loop almost unchanged: a model with parameters θ\theta makes predictions, a loss compares them with labels, and an optimiser moves θ\theta against the gradient of the average loss. In Stage 4 the model will be a parameterised quantum circuit and the gradient will be estimated from measurement shots; in this lesson every object is classical, so each step of the loop can be derived exactly. We work with real vectors x∈Rdx\in\mathbb R^d, the Euclidean inner product u⊤v=∑jujvju^\top v=\sum_ju_jv_j and norm ∥u∥=u⊤u\lVert u\rVert=\sqrt{u^\top u} — the real special case of the Hilbert-space structure of the first lessons.

Definition

Supervised learning problem

An unknown probability distribution D\mathcal D on X×Y\mathcal X\times\mathcal Y generates labelled examples. A training sample S=((x1,y1),…,(xm,ym))S=\big((x_1,y_1),\dots,(x_m,y_m)\big) consists of mm independent draws from D\mathcal D. A hypothesis class F\mathcal F is a set of functions f:X→Y^f:\mathcal X\to\hat{\mathcal Y}, and a loss ℓ:Y^×Y→[0,∞)\ell:\hat{\mathcal Y}\times\mathcal Y\to[0,\infty) scores a prediction against a label: the squared loss ℓ(y^,y)=(y^−y)2\ell(\hat y,y)=(\hat y-y)^2 for regression, the cross-entropy ℓ(p,y)=−ylog⁡p−(1−y)log⁡(1−p)\ell(p,y)=-y\log p-(1-y)\log(1-p) for binary classification with y∈{0,1}y\in\{0,1\}. The population risk of ff is R(f)=E(x,y)∼D ℓ(f(x),y)R(f)=\mathbb E_{(x,y)\sim\mathcal D}\,\ell(f(x),y) and its empirical risk on SS is R^S(f)=1m∑i=1mℓ(f(xi),yi)\hat R_S(f)=\tfrac1m\sum_{i=1}^m\ell(f(x_i),y_i).

Definition

Empirical-risk minimisation and generalisation gap

Empirical-risk minimisation returns a hypothesis with the smallest training loss, f^S∈arg min⁡f∈F 1m∑i=1mℓ(f(xi),yi).\hat f_S\in\operatorname*{arg\,min}_{f\in\mathcal F}\ \frac1m\sum_{i=1}^m\ell\big(f(x_i),y_i\big). The population risk cannot be computed because D\mathcal D is unknown; ERM replaces it by its sample average. The generalisation gap of the returned hypothesis is R(f^S)−R^S(f^S)R(\hat f_S)-\hat R_S(\hat f_S), the amount by which the training loss under-reports the loss on fresh data. When F={fθ}\mathcal F=\{f_\theta\} is parameterised, ERM is the finite-dimensional optimisation problem min⁡θL(θ)\min_\theta\mathcal L(\theta) with L(θ)=R^S(fθ)\mathcal L(\theta)=\hat R_S(f_\theta), which the rest of the lesson solves by gradient descent.

Proposition

Training risk is unbiased for a fixed model and optimistic for the ERM model

Assume all risks below are finite. (a) For every fixed ff that does not depend on SS, ES[R^S(f)]=R(f)\mathbb E_S[\hat R_S(f)]=R(f). (b) If f^S\hat f_S is an empirical-risk minimiser over F\mathcal F, then ES[R^S(f^S)]≤inf⁡f∈FR(f)≤ES[R(f^S)],\mathbb E_S\big[\hat R_S(\hat f_S)\big]\le\inf_{f\in\mathcal F}R(f)\le\mathbb E_S\big[R(\hat f_S)\big], so the expected generalisation gap ES[R(f^S)−R^S(f^S)]\mathbb E_S[R(\hat f_S)-\hat R_S(\hat f_S)] is non-negative.

Proof (7 steps)
  1. Step 1:
    ES[R^S(f)]=1m∑i=1mE[ℓ(f(xi),yi)]\mathbb E_S\big[\hat R_S(f)\big]=\frac1m\sum_{i=1}^m\mathbb E\big[\ell(f(x_i),y_i)\big]
    Linearity of expectation applied to the average of mm terms.
  2. Step 2:
    E[ℓ(f(xi),yi)]=R(f)  ∀i⟹ES[R^S(f)]=1m⋅mR(f)=R(f)\mathbb E\big[\ell(f(x_i),y_i)\big]=R(f)\ \ \forall i\quad\Longrightarrow\quad\mathbb E_S\big[\hat R_S(f)\big]=\frac1m\cdot mR(f)=R(f)
    Each pair (xi,yi)(x_i,y_i) has distribution D\mathcal D, and ff is fixed, so every term has the expectation R(f)R(f). This proves (a).
  3. Step 3:
    R^S(f^S)≤R^S(f)∀f∈F, ∀S\hat R_S(\hat f_S)\le\hat R_S(f)\qquad\forall f\in\mathcal F,\ \forall S
    Definition of f^S\hat f_S as a minimiser of R^S\hat R_S over F\mathcal F.
  4. Step 4:
    ES[R^S(f^S)]≤ES[R^S(f)]=R(f)∀f∈F\mathbb E_S\big[\hat R_S(\hat f_S)\big]\le\mathbb E_S\big[\hat R_S(f)\big]=R(f)\qquad\forall f\in\mathcal F
    Expectation preserves inequalities; then apply (a) to the fixed ff.
  5. Step 5:
    ES[R^S(f^S)]≤inf⁡f∈FR(f)\mathbb E_S\big[\hat R_S(\hat f_S)\big]\le\inf_{f\in\mathcal F}R(f)
    The left side does not depend on ff, so it is a lower bound of the set {R(f):f∈F}\{R(f):f\in\mathcal F\}, hence below its infimum.
  6. Step 6:
    R(f^S)≥inf⁡f∈FR(f)  ∀S⟹ES[R(f^S)]≥inf⁡f∈FR(f)R(\hat f_S)\ge\inf_{f\in\mathcal F}R(f)\ \ \forall S\quad\Longrightarrow\quad\mathbb E_S\big[R(\hat f_S)\big]\ge\inf_{f\in\mathcal F}R(f)
    f^S∈F\hat f_S\in\mathcal F for every sample, and the expectation of a quantity bounded below by a constant is bounded below by that constant.
  7. Step 7:
    ES[R(f^S)−R^S(f^S)]≥inf⁡fR(f)−inf⁡fR(f)=0\mathbb E_S\big[R(\hat f_S)-\hat R_S(\hat f_S)\big]\ge\inf_{f}R(f)-\inf_fR(f)=0
    Subtract the two previous inequalities. Part (a) fails for f^S\hat f_S precisely because f^S\hat f_S was chosen by looking at SS.

∎

For regression with the squared loss, the error on a fresh example can be split further. Fix a test input xx and write the label as y=f∗(x)+εy=f^\ast(x)+\varepsilon with the regression function f∗(x)=E[y∣x]f^\ast(x)=\mathbb E[y\mid x] and noise ε\varepsilon satisfying E[ε∣x]=0\mathbb E[\varepsilon\mid x]=0 and E[ε2∣x]=σ2(x)\mathbb E[\varepsilon^2\mid x]=\sigma^2(x). The prediction f^S(x)\hat f_S(x) is random through the training sample SS, which is independent of the test noise ε\varepsilon. Let fˉ(x)=ES[f^S(x)]\bar f(x)=\mathbb E_S[\hat f_S(x)] be the average prediction over training samples.

Theorem

Bias–variance decomposition of the expected squared error

Under the assumptions above, for every input xx, ES,ε[(y−f^S(x))2]=σ2(x)⏟noise+(f∗(x)−fˉ(x))2⏟bias2+ES[(f^S(x)−fˉ(x))2]⏟variance,\mathbb E_{S,\varepsilon}\big[(y-\hat f_S(x))^2\big]=\underbrace{\sigma^2(x)}_{\text{noise}}+\underbrace{\big(f^\ast(x)-\bar f(x)\big)^2}_{\text{bias}^2}+\underbrace{\mathbb E_S\big[(\hat f_S(x)-\bar f(x))^2\big]}_{\text{variance}}, and averaging over x∼Dxx\sim\mathcal D_x gives the same decomposition of the expected test MSE.

Proof (7 steps)
  1. Step 1:
    y−f^S(x)=ε+(f∗(x)−fˉ(x))⏟=:a+(fˉ(x)−f^S(x))⏟=:By-\hat f_S(x)=\varepsilon+\underbrace{\big(f^\ast(x)-\bar f(x)\big)}_{=:a}+\underbrace{\big(\bar f(x)-\hat f_S(x)\big)}_{=:B}
    Substitute y=f∗(x)+εy=f^\ast(x)+\varepsilon, then add and subtract fˉ(x)=ES[f^S(x)]\bar f(x)=\mathbb E_S[\hat f_S(x)]. At fixed xx, aa is a number and BB is random only through SS.
  2. Step 2:
    (y−f^S(x))2=ε2+a2+B2+2aε+2εB+2aB\big(y-\hat f_S(x)\big)^2=\varepsilon^2+a^2+B^2+2a\varepsilon+2\varepsilon B+2aB
    Expand (u+v+w)2=u2+v2+w2+2uv+2uw+2vw(u+v+w)^2=u^2+v^2+w^2+2uv+2uw+2vw with u=εu=\varepsilon, v=av=a, w=Bw=B.
  3. Step 3:
    E[2aε]=2a E[ε∣x]=0\mathbb E[2a\varepsilon]=2a\,\mathbb E[\varepsilon\mid x]=0
    aa is constant at fixed xx and the noise has conditional mean zero.
  4. Step 4:
    E[2εB]=2 E[ε∣x]  ES[B]=0\mathbb E[2\varepsilon B]=2\,\mathbb E[\varepsilon\mid x]\;\mathbb E_S[B]=0
    The test noise is independent of the training sample, so the expectation of the product factorises; the first factor is 00.
  5. Step 5:
    E[2aB]=2a(fˉ(x)−ES[f^S(x)])=2a⋅0=0\mathbb E[2aB]=2a\big(\bar f(x)-\mathbb E_S[\hat f_S(x)]\big)=2a\cdot0=0
    Linearity of ES\mathbb E_S and the definition of fˉ(x)\bar f(x): the cross term vanishes because BB has mean zero.
  6. Step 6:
    ES,ε[(y−f^S(x))2]=E[ε2∣x]+a2+ES[B2]=σ2(x)+(f∗(x)−fˉ(x))2+ES[(f^S(x)−fˉ(x))2]\mathbb E_{S,\varepsilon}\big[(y-\hat f_S(x))^2\big]=\mathbb E[\varepsilon^2\mid x]+a^2+\mathbb E_S[B^2]=\sigma^2(x)+\big(f^\ast(x)-\bar f(x)\big)^2+\mathbb E_S\big[(\hat f_S(x)-\bar f(x))^2\big]
    Take expectations of the expansion, drop the three vanishing cross terms, and use B2=(f^S(x)−fˉ(x))2B^2=(\hat f_S(x)-\bar f(x))^2.
  7. Step 7:
    Ex,S,ε[(y−f^S(x))2]=Ex[σ2(x)]+Ex[(f∗(x)−fˉ(x))2]+Ex ES[(f^S(x)−fˉ(x))2]\mathbb E_{x,S,\varepsilon}\big[(y-\hat f_S(x))^2\big]=\mathbb E_x[\sigma^2(x)]+\mathbb E_x\big[(f^\ast(x)-\bar f(x))^2\big]+\mathbb E_x\,\mathbb E_S\big[(\hat f_S(x)-\bar f(x))^2\big]
    Tower property: average the pointwise identity over the test input x∼Dxx\sim\mathcal D_x, which is independent of SS.

∎

Quick check

A regression model reaches training MSE R^S(f^S)=0.01\hat R_S(\hat f_S)=0.01 but test MSE 0.400.40 on a large held-out set drawn from the same distribution. Which statement is supported by this lesson?

Select one answer

Linear regression and the gradient of the mean-squared error#

The linear model predicts y^=w⊤x+b\hat y=w^\top x+b with weights w∈Rdw\in\mathbb R^d and bias b∈Rb\in\mathbb R. Collect the inputs as the rows of the design matrix X∈Rm×dX\in\mathbb R^{m\times d} (so XikX_{ik} is feature kk of sample ii), the labels as y∈Rmy\in\mathbb R^m, and write 1∈Rm\mathbf 1\in\mathbb R^m for the all-ones vector. ERM with the squared loss minimises the mean-squared error L(w,b)=1m∑i=1m(y^i−yi)2,y^i=w⊤xi+b=∑j=1dwjXij+b,\mathcal L(w,b)=\frac1m\sum_{i=1}^m\big(\hat y_i-y_i\big)^2,\qquad \hat y_i=w^\top x_i+b=\sum_{j=1}^dw_jX_{ij}+b, and the residuals ri=y^i−yir_i=\hat y_i-y_i form the vector r=Xw+b1−yr=Xw+b\mathbf 1-y.

Derivation

Gradient of the MSE, coordinate by coordinate

  1. Step 1:
    ∂L∂wk=1m∑i=1m∂∂wk ri2\frac{\partial\mathcal L}{\partial w_k}=\frac1m\sum_{i=1}^m\frac{\partial}{\partial w_k}\,r_i^2
    Differentiation is linear, so it passes through the constant factor 1m\tfrac1m and the finite sum.
  2. Step 2:
    ∂∂wk ri2=2ri ∂ri∂wk\frac{\partial}{\partial w_k}\,r_i^2=2r_i\,\frac{\partial r_i}{\partial w_k}
    Chain rule with the outer function u↦u2u\mapsto u^2, whose derivative is 2u2u.
  3. Step 3:
    ∂ri∂wk=∂∂wk(∑j=1dwjXij+b−yi)=Xik\frac{\partial r_i}{\partial w_k}=\frac{\partial}{\partial w_k}\Big(\sum_{j=1}^dw_jX_{ij}+b-y_i\Big)=X_{ik}
    Only the term j=kj=k of the sum depends on wkw_k; bb and yiy_i are constants with respect to wkw_k.
  4. Step 4:
    ∂L∂wk=2m∑i=1mri Xik,k=1,…,d\frac{\partial\mathcal L}{\partial w_k}=\frac2m\sum_{i=1}^m r_i\,X_{ik},\qquad k=1,\dots,d
    Substitute the two previous lines into the first.
  5. Step 5:
    ∇wL=2m∑i=1m(y^i−yi) xi\nabla_w\mathcal L=\frac2m\sum_{i=1}^m\big(\hat y_i-y_i\big)\,x_i
    Stack the dd coordinates: the vector with entries XikX_{ik}, k=1,…,dk=1,\dots,d, is xix_i.
  6. Step 6:
    ∂ri∂b=1⟹∂L∂b=2m∑i=1m(y^i−yi)\frac{\partial r_i}{\partial b}=1\quad\Longrightarrow\quad\frac{\partial\mathcal L}{\partial b}=\frac2m\sum_{i=1}^m\big(\hat y_i-y_i\big)
    Same chain rule with ∂ri/∂b=1\partial r_i/\partial b=1, because bb enters every prediction with coefficient 11.

Result

∇wL=2m∑i=1m(y^i−yi) xi,∂L∂b=2m∑i=1m(y^i−yi)\boxed{\nabla_w\mathcal L=\frac2m\sum_{i=1}^m(\hat y_i-y_i)\,x_i,\qquad \frac{\partial\mathcal L}{\partial b}=\frac2m\sum_{i=1}^m(\hat y_i-y_i)}
Derivation

Matrix form and a worked example

  1. Step 1:
    (X⊤r)k=∑i=1m(X⊤)ki ri=∑i=1mXik ri\big(X^\top r\big)_k=\sum_{i=1}^m(X^\top)_{ki}\,r_i=\sum_{i=1}^mX_{ik}\,r_i
    Definition of the matrix–vector product and of the transpose, (X⊤)ki=Xik(X^\top)_{ki}=X_{ik}.
  2. Step 2:
    ∇wL=2mX⊤r,∂L∂b=2m1⊤r,r=Xw+b1−y\nabla_w\mathcal L=\frac2mX^\top r,\qquad \frac{\partial\mathcal L}{\partial b}=\frac2m\mathbf 1^\top r,\qquad r=Xw+b\mathbf 1-y
    Compare with the coordinate formula: entry kk of X⊤rX^\top r is the dot product of column kk of XX with rr, and 1⊤r=∑iri\mathbf 1^\top r=\sum_ir_i.
  3. Step 3:
    w=(0.5−0.25), b=0.1, X=(1234), y=(12):Xw+b1=(0.5−0.5+0.11.5−1+0.1)=(0.10.6)w=\begin{pmatrix}0.5\\-0.25\end{pmatrix},\ b=0.1,\ X=\begin{pmatrix}1&2\\3&4\end{pmatrix},\ y=\begin{pmatrix}1\\2\end{pmatrix}:\quad Xw+b\mathbf 1=\begin{pmatrix}0.5-0.5+0.1\\1.5-1+0.1\end{pmatrix}=\begin{pmatrix}0.1\\0.6\end{pmatrix}
    The predictions of the example (m=d=2m=d=2) used as the visible test of the lab.
  4. Step 4:
    r=(0.1−10.6−2)=(−0.9−1.4)r=\begin{pmatrix}0.1-1\\0.6-2\end{pmatrix}=\begin{pmatrix}-0.9\\-1.4\end{pmatrix}
    Residuals are prediction minus label.
  5. Step 5:
    ∇wL=22(1324)(−0.9−1.4)=(−0.9−4.2−1.8−5.6)=(−5.1−7.4)\nabla_w\mathcal L=\frac22\begin{pmatrix}1&3\\2&4\end{pmatrix}\begin{pmatrix}-0.9\\-1.4\end{pmatrix}=\begin{pmatrix}-0.9-4.2\\-1.8-5.6\end{pmatrix}=\begin{pmatrix}-5.1\\-7.4\end{pmatrix}
    Apply 2mX⊤r\tfrac2mX^\top r with m=2m=2, so the prefactor is 11.
  6. Step 6:
    ∂L∂b=22 (−0.9−1.4)=−2.3\frac{\partial\mathcal L}{\partial b}=\frac22\,(-0.9-1.4)=-2.3
    Apply 2m1⊤r\tfrac2m\mathbf 1^\top r. All components are negative, so increasing w1w_1, w2w_2 and bb lowers the loss: both predictions are too small.

Result

∇wL=2mX⊤(Xw+b1−y),∂bL=2m1⊤(Xw+b1−y)\boxed{\nabla_w\mathcal L=\tfrac2mX^\top(Xw+b\mathbf 1-y),\qquad \partial_b\mathcal L=\tfrac2m\mathbf 1^\top(Xw+b\mathbf 1-y)}
Derivation

The MSE is convex and L-smooth

  1. Step 1:
    X~=[ X  1 ]∈Rm×(d+1),θ=(wb)⟹r=X~θ−y,∇θL=2mX~⊤(X~θ−y)\tilde X=\big[\,X\ \ \mathbf 1\,\big]\in\mathbb R^{m\times(d+1)},\quad \theta=\begin{pmatrix}w\\b\end{pmatrix}\quad\Longrightarrow\quad r=\tilde X\theta-y,\quad \nabla_\theta\mathcal L=\frac2m\tilde X^\top(\tilde X\theta-y)
    Append a constant feature 11 so that the bias becomes an ordinary weight; the previous derivation then covers ww and bb at once.
  2. Step 2:
    ∂2L∂θk ∂θl=2m∑i=1mX~ikX~il⟹∇2L=2mX~⊤X~\frac{\partial^2\mathcal L}{\partial\theta_k\,\partial\theta_l}=\frac2m\sum_{i=1}^m\tilde X_{ik}\tilde X_{il}\quad\Longrightarrow\quad\nabla^2\mathcal L=\frac2m\tilde X^\top\tilde X
    Differentiate 2m∑iX~ikri\tfrac2m\sum_i\tilde X_{ik}r_i with respect to θl\theta_l using ∂ri/∂θl=X~il\partial r_i/\partial\theta_l=\tilde X_{il}. The Hessian does not depend on θ\theta.
  3. Step 3:
    v⊤∇2L v=2m∥X~v∥2≥0∀v∈Rd+1v^\top\nabla^2\mathcal L\,v=\frac2m\lVert\tilde Xv\rVert^2\ge0\qquad\forall v\in\mathbb R^{d+1}
    v⊤X~⊤X~v=(X~v)⊤(X~v)v^\top\tilde X^\top\tilde Xv=(\tilde Xv)^\top(\tilde Xv). A twice-differentiable function with positive-semidefinite Hessian everywhere is convex.
  4. Step 4:
    ∥∇L(θ)−∇L(θ′)∥=∥2mX~⊤X~(θ−θ′)∥≤λmax⁡(2mX~⊤X~)∥θ−θ′∥\lVert\nabla\mathcal L(\theta)-\nabla\mathcal L(\theta')\rVert=\Big\lVert\frac2m\tilde X^\top\tilde X(\theta-\theta')\Big\rVert\le\lambda_{\max}\Big(\frac2m\tilde X^\top\tilde X\Big)\lVert\theta-\theta'\rVert
    The gradient is affine in θ\theta, and a symmetric positive-semidefinite matrix stretches no vector by more than its largest eigenvalue (spectral theorem).

Result

L is convex and L-smooth with L=2m λmax⁡(X~⊤X~)\boxed{\mathcal L\ \text{is convex and}\ L\text{-smooth with}\ L=\frac2m\,\lambda_{\max}\big(\tilde X^\top\tilde X\big)}

Quick check

Keep w=(0.5,−0.25)w=(0.5,-0.25), X=(1234)X=\begin{pmatrix}1&2\\3&4\end{pmatrix}, y=(1,2)y=(1,2) from the worked example but set b=1.25b=1.25. What is ∂L/∂b\partial\mathcal L/\partial b?

tolerance ±1e-9

The listing below implements the vectorised gradient, checks it against central finite differences (L(θ+hej)−L(θ−hej))/(2h)\big(\mathcal L(\theta+he_j)-\mathcal L(\theta-he_j)\big)/(2h), computes L=2mλmax⁡(X~⊤X~)L=\tfrac2m\lambda_{\max}(\tilde X^\top\tilde X) by power iteration and runs gradient descent with step 1/L1/L on a noiseless problem whose minimiser is w∗=(2,−1)w^\ast=(2,-1), b∗=0.5b^\ast=0.5 with L∗=0\mathcal L^\ast=0. Each printed loss is compared with the bound L∥θ0−θ∗∥2/(2k)L\lVert\theta_0-\theta^\ast\rVert^2/(2k) proved in the convergence section, with ∥θ0−θ∗∥2=4+1+0.25=5.25\lVert\theta_0-\theta^\ast\rVert^2=4+1+0.25=5.25.

TypeScriptrunnable
// Linear least squares with two features: yhat = w . x + b, loss = (1/m) sum_i (yhat_i - y_i)^2.
const X = [[1, 2], [2, 0], [3, 1], [0, 1], [-1, 3], [2, 2]];
const y = X.map(([x1, x2]) => 2 * x1 - x2 + 0.5); // noiseless: w* = (2, -1), b* = 0.5, loss* = 0
const m = X.length;
const residuals = (w, b) => X.map((x, i) => x[0] * w[0] + x[1] * w[1] + b - y[i]);
const loss = (w, b) => residuals(w, b).reduce((s, r) => s + r * r, 0) / m;

// Vectorised gradient: grad_w = (2/m) X^T r and dL/db = (2/m) 1^T r with r = Xw + b1 - y.
function gradient(w, b) {
  const r = residuals(w, b);
  const gw = [0, 1].map((k) => (2 / m) * X.reduce((s, x, i) => s + x[k] * r[i], 0));
  const gb = (2 / m) * r.reduce((s, ri) => s + ri, 0);
  return { gw, gb };
}

// Check the gradient against central finite differences at an arbitrary point.
const w0 = [0.3, -0.7], b0 = 0.2, h = 1e-6;
const g0 = gradient(w0, b0);
const fd = [
  (loss([w0[0] + h, w0[1]], b0) - loss([w0[0] - h, w0[1]], b0)) / (2 * h),
  (loss([w0[0], w0[1] + h], b0) - loss([w0[0], w0[1] - h], b0)) / (2 * h),
  (loss(w0, b0 + h) - loss(w0, b0 - h)) / (2 * h),
];
const fdErr = Math.max(...[...g0.gw, g0.gb].map((g, j) => Math.abs(g - fd[j])));
console.log('finite-difference check: ' + (fdErr < 1e-6 ? 'ok' : 'FAILED'));

// Smoothness constant L = lambda_max of the Hessian (2/m) A^T A with A = [X 1], by power iteration.
const A = X.map((x) => [x[0], x[1], 1]);
const hess = [0, 1, 2].map((p) => [0, 1, 2].map((q) => (2 / m) * A.reduce((s, a) => s + a[p] * a[q], 0)));
let v = [1, 1, 1];
for (let it = 0; it < 500; it++) {
  const hv = hess.map((row) => row.reduce((s, hpq, q) => s + hpq * v[q], 0));
  const nrm = Math.hypot(...hv);
  v = hv.map((t) => t / nrm);
}
const L = hess.map((row) => row.reduce((s, hpq, q) => s + hpq * v[q], 0)).reduce((s, t, p) => s + t * v[p], 0);
console.log('L = ' + L.toFixed(6));

// Gradient descent with step 1/L from (w, b) = (0, 0, 0); bound L * ||x0 - x*||^2 / (2k).
let w = [0, 0], b = 0;
const dist0 = 2 ** 2 + (-1) ** 2 + 0.5 ** 2;
const report = [1, 2, 5, 10, 20, 50, 100, 200, 400];
console.log('k = 0  loss = ' + loss(w, b).toExponential(4));
for (let k = 1; k <= 400; k++) {
  const { gw, gb } = gradient(w, b);
  w = [w[0] - gw[0] / L, w[1] - gw[1] / L];
  b = b - gb / L;
  if (report.includes(k)) {
    const bound = (L * dist0) / (2 * k);
    console.log('k = ' + k + '  loss = ' + loss(w, b).toExponential(4) + '  bound = ' + bound.toExponential(4) + (loss(w, b) <= bound ? '  ok' : '  VIOLATED'));
  }
}
console.log('w = (' + w[0].toFixed(6) + ', ' + w[1].toFixed(6) + ')  b = ' + b.toFixed(6));

Expected output

finite-difference check: ok
L = 10.097998
k = 0  loss = 1.2917e+1
k = 1  loss = 3.1800e+0  bound = 2.6507e+1  ok
k = 2  loss = 1.0315e+0  bound = 1.3254e+1  ok
k = 5  loss = 4.5253e-2  bound = 5.3014e+0  ok
k = 10  loss = 9.4194e-3  bound = 2.6507e+0  ok
k = 20  loss = 6.2597e-3  bound = 1.3254e+0  ok
k = 50  loss = 1.9073e-3  bound = 5.3014e-1  ok
k = 100  loss = 2.6314e-4  bound = 2.6507e-1  ok
k = 200  loss = 5.0089e-6  bound = 1.3254e-1  ok
k = 400  loss = 1.8148e-9  bound = 6.6268e-2  ok
w = (2.000030, -0.999950)  b = 0.499878

Logistic regression and the cross-entropy gradient#

For labels yi∈{0,1}y_i\in\{0,1\}, logistic regression models the probability of class 11 as pi=σ(zi)p_i=\sigma(z_i) with the logit zi=w⊤xi+bz_i=w^\top x_i+b and the sigmoid σ(z)=11+e−z∈(0,1)\sigma(z)=\dfrac1{1+e^{-z}}\in(0,1). Maximising the likelihood ∏ipiyi(1−pi)1−yi\prod_ip_i^{y_i}(1-p_i)^{1-y_i} is the same as minimising the average negative log-likelihood, the binary cross-entropy L(w,b)=1m∑i=1mℓi,ℓi=−yilog⁡pi−(1−yi)log⁡(1−pi).\mathcal L(w,b)=\frac1m\sum_{i=1}^m\ell_i,\qquad \ell_i=-y_i\log p_i-(1-y_i)\log(1-p_i).

Derivation

Derivative of the sigmoid

  1. Step 1:
    σ(z)=(1+e−z)−1\sigma(z)=\big(1+e^{-z}\big)^{-1}
    Definition, written as a power so that the chain rule applies directly.
  2. Step 2:
    σ′(z)=−(1+e−z)−2⋅(−e−z)=e−z(1+e−z)2\sigma'(z)=-\big(1+e^{-z}\big)^{-2}\cdot\big(-e^{-z}\big)=\frac{e^{-z}}{(1+e^{-z})^2}
    Chain rule: the derivative of u−1u^{-1} is −u−2-u^{-2}, and ddz(1+e−z)=−e−z\tfrac{d}{dz}(1+e^{-z})=-e^{-z}.
  3. Step 3:
    e−z(1+e−z)2=11+e−z⋅e−z1+e−z\frac{e^{-z}}{(1+e^{-z})^2}=\frac1{1+e^{-z}}\cdot\frac{e^{-z}}{1+e^{-z}}
    Split the square in the denominator into two factors.
  4. Step 4:
    e−z1+e−z=(1+e−z)−11+e−z=1−11+e−z=1−σ(z)\frac{e^{-z}}{1+e^{-z}}=\frac{(1+e^{-z})-1}{1+e^{-z}}=1-\frac1{1+e^{-z}}=1-\sigma(z)
    Add and subtract 11 in the numerator.
  5. Step 5:
    σ′(z)=σ(z)(1−σ(z))\sigma'(z)=\sigma(z)\big(1-\sigma(z)\big)
    Combine the two previous lines. In particular σ′(0)=12⋅12=14\sigma'(0)=\tfrac12\cdot\tfrac12=\tfrac14, the maximum of σ′\sigma', since p(1−p)≤14p(1-p)\le\tfrac14 with equality at p=12p=\tfrac12.

Result

σ′(z)=σ(z)(1−σ(z))≤14\boxed{\sigma'(z)=\sigma(z)\big(1-\sigma(z)\big)\le\tfrac14}
Derivation

Gradient of the binary cross-entropy

  1. Step 1:
    ∂ℓi∂pi=−yipi+1−yi1−pi\frac{\partial\ell_i}{\partial p_i}=-\frac{y_i}{p_i}+\frac{1-y_i}{1-p_i}
    ddplog⁡p=1p\tfrac{d}{dp}\log p=\tfrac1p and ddplog⁡(1−p)=−11−p\tfrac{d}{dp}\log(1-p)=-\tfrac1{1-p}.
  2. Step 2:
    ∂ℓi∂pi=−yi(1−pi)+(1−yi)pipi(1−pi)=pi−yipi(1−pi)\frac{\partial\ell_i}{\partial p_i}=\frac{-y_i(1-p_i)+(1-y_i)p_i}{p_i(1-p_i)}=\frac{p_i-y_i}{p_i(1-p_i)}
    Common denominator pi(1−pi)p_i(1-p_i); the numerator simplifies as −yi+yipi+pi−yipi=pi−yi-y_i+y_ip_i+p_i-y_ip_i=p_i-y_i.
  3. Step 3:
    ∂pi∂zi=σ′(zi)=pi(1−pi)\frac{\partial p_i}{\partial z_i}=\sigma'(z_i)=p_i(1-p_i)
    The sigmoid derivative derived above, with pi=σ(zi)p_i=\sigma(z_i).
  4. Step 4:
    ∂ℓi∂zi=∂ℓi∂pi ∂pi∂zi=pi−yipi(1−pi)⋅pi(1−pi)=σ(zi)−yi\frac{\partial\ell_i}{\partial z_i}=\frac{\partial\ell_i}{\partial p_i}\,\frac{\partial p_i}{\partial z_i}=\frac{p_i-y_i}{p_i(1-p_i)}\cdot p_i(1-p_i)=\sigma(z_i)-y_i
    Chain rule; the factor pi(1−pi)∈(0,14]p_i(1-p_i)\in(0,\tfrac14] is non-zero and cancels.
  5. Step 5:
    ∂zi∂wk=Xik,∂zi∂b=1\frac{\partial z_i}{\partial w_k}=X_{ik},\qquad\frac{\partial z_i}{\partial b}=1
    The logit zi=∑jwjXij+bz_i=\sum_jw_jX_{ij}+b is linear in the parameters, exactly as the prediction of the linear model.
  6. Step 6:
    ∇wL=1m∑i=1m(σ(zi)−yi) xi=1mX⊤(p−y),∂L∂b=1m∑i=1m(σ(zi)−yi)\nabla_w\mathcal L=\frac1m\sum_{i=1}^m\big(\sigma(z_i)-y_i\big)\,x_i=\frac1mX^\top(p-y),\qquad\frac{\partial\mathcal L}{\partial b}=\frac1m\sum_{i=1}^m\big(\sigma(z_i)-y_i\big)
    Average the per-sample chain rule ∂ℓi/∂wk=(σ(zi)−yi)Xik\partial\ell_i/\partial w_k=(\sigma(z_i)-y_i)X_{ik} over ii and stack the coordinates.

Result

∇wL=1m∑i=1m(σ(zi)−yi) xi\boxed{\nabla_w\mathcal L=\frac1m\sum_{i=1}^m\big(\sigma(z_i)-y_i\big)\,x_i}
Derivation

The cross-entropy is convex and L-smooth

  1. Step 1:
    ∂2L∂θl ∂θk=1m∑i=1m∂σ(zi)∂θl X~ik=1m∑i=1mpi(1−pi) X~ilX~ik\frac{\partial^2\mathcal L}{\partial\theta_l\,\partial\theta_k}=\frac1m\sum_{i=1}^m\frac{\partial\sigma(z_i)}{\partial\theta_l}\,\tilde X_{ik}=\frac1m\sum_{i=1}^mp_i(1-p_i)\,\tilde X_{il}\tilde X_{ik}
    Differentiate ∂θkL=1m∑i(σ(zi)−yi)X~ik\partial_{\theta_k}\mathcal L=\tfrac1m\sum_i(\sigma(z_i)-y_i)\tilde X_{ik} with θ=(w,b)\theta=(w,b) and X~=[X 1]\tilde X=[X\ \mathbf 1] (rows x~i⊤=(xi⊤,1)\tilde x_i^\top=(x_i^\top,1)), and use ∂σ(zi)/∂θl=pi(1−pi)X~il\partial\sigma(z_i)/\partial\theta_l=p_i(1-p_i)\tilde X_{il}.
  2. Step 2:
    v⊤∇2L v=1m∑i=1mpi(1−pi)(x~i⊤v)2≥0v^\top\nabla^2\mathcal L\,v=\frac1m\sum_{i=1}^mp_i(1-p_i)\big(\tilde x_i^\top v\big)^2\ge0
    Each weight pi(1−pi)p_i(1-p_i) lies in (0,14](0,\tfrac14] and each square is non-negative, so the Hessian is positive semidefinite and L\mathcal L is convex.
  3. Step 3:
    v⊤∇2L v≤14m∑i=1m(x~i⊤v)2=14m∥X~v∥2≤λmax⁡(X~⊤X~)4m∥v∥2v^\top\nabla^2\mathcal L\,v\le\frac1{4m}\sum_{i=1}^m\big(\tilde x_i^\top v\big)^2=\frac1{4m}\lVert\tilde Xv\rVert^2\le\frac{\lambda_{\max}(\tilde X^\top\tilde X)}{4m}\lVert v\rVert^2
    Bound pi(1−pi)≤14p_i(1-p_i)\le\tfrac14, then use the Rayleigh quotient bound of the spectral theorem.

Result

LBCE is convex and L-smooth with L≤λmax⁡(X~⊤X~)4m\boxed{\mathcal L_{\mathrm{BCE}}\ \text{is convex and}\ L\text{-smooth with}\ L\le\frac{\lambda_{\max}(\tilde X^\top\tilde X)}{4m}}

Quick check

At a sample with logit zi=0z_i=0 and label yi=1y_i=1, what is the per-sample derivative ∂ℓi/∂zi\partial\ell_i/\partial z_i of the binary cross-entropy?

tolerance ±1e-9

Smoothness and the descent lemma#

Gradient descent (GD) with step size (learning rate) η>0\eta>0 iterates xk+1=xk−η ∇f(xk).x_{k+1}=x_k-\eta\,\nabla f(x_k). The first-order Taylor model f(xk)−η∥∇f(xk)∥2f(x_k)-\eta\lVert\nabla f(x_k)\rVert^2 predicts a decrease, but only for small enough η\eta; to control the error of the linear model we need a bound on how fast the gradient can change.

Definition

L-smooth function

A differentiable function f:Rn→Rf:\mathbb R^n\to\mathbb R is LL-smooth if its gradient is LL-Lipschitz: ∥∇f(x)−∇f(y)∥≤L ∥x−y∥∀x,y∈Rn.\lVert\nabla f(x)-\nabla f(y)\rVert\le L\,\lVert x-y\rVert\qquad\forall x,y\in\mathbb R^n. Example: for f(x)=12x⊤Axf(x)=\tfrac12x^\top Ax with AA symmetric positive semidefinite, ∇f(x)=Ax\nabla f(x)=Ax and ∥A(x−y)∥≤λmax⁡(A)∥x−y∥\lVert A(x-y)\rVert\le\lambda_{\max}(A)\lVert x-y\rVert, so L=λmax⁡(A)L=\lambda_{\max}(A). The MSE and the cross-entropy are LL-smooth with the constants derived above.

Lemma

Descent lemma

If ff is LL-smooth, then for all x,y∈Rnx,y\in\mathbb R^n f(y)≤f(x)+∇f(x)⊤(y−x)+L2∥y−x∥2.f(y)\le f(x)+\nabla f(x)^\top(y-x)+\frac L2\lVert y-x\rVert^2 .

Proof (6 steps)
  1. Step 1:
    φ(t)=f(x+t(y−x)),φ′(t)=∇f(x+t(y−x))⊤(y−x)\varphi(t)=f\big(x+t(y-x)\big),\qquad\varphi'(t)=\nabla f\big(x+t(y-x)\big)^\top(y-x)
    Restrict ff to the segment from xx to yy; the multivariate chain rule gives φ′\varphi'. A Lipschitz gradient is continuous, so φ\varphi is continuously differentiable.
  2. Step 2:
    f(y)−f(x)=φ(1)−φ(0)=∫01∇f(x+t(y−x))⊤(y−x) dtf(y)-f(x)=\varphi(1)-\varphi(0)=\int_0^1\nabla f\big(x+t(y-x)\big)^\top(y-x)\,dt
    Fundamental theorem of calculus for φ\varphi on [0,1][0,1] (the integral form of Taylor's theorem of order one).
  3. Step 3:
    f(y)−f(x)−∇f(x)⊤(y−x)=∫01(∇f(x+t(y−x))−∇f(x))⊤(y−x) dtf(y)-f(x)-\nabla f(x)^\top(y-x)=\int_0^1\Big(\nabla f\big(x+t(y-x)\big)-\nabla f(x)\Big)^\top(y-x)\,dt
    Subtract ∇f(x)⊤(y−x)=∫01∇f(x)⊤(y−x) dt\nabla f(x)^\top(y-x)=\int_0^1\nabla f(x)^\top(y-x)\,dt, the integral of a constant over an interval of length 11.
  4. Step 4:
    ≤∫01∥∇f(x+t(y−x))−∇f(x)∥ ∥y−x∥ dt\le\int_0^1\Big\lVert\nabla f\big(x+t(y-x)\big)-\nabla f(x)\Big\rVert\,\lVert y-x\rVert\,dt
    Cauchy–Schwarz u⊤v≤∥u∥∥v∥u^\top v\le\lVert u\rVert\lVert v\rVert inside the integral; integrals preserve pointwise inequalities.
  5. Step 5:
    ≤∫01L ∥t(y−x)∥ ∥y−x∥ dt=L∥y−x∥2∫01t dt\le\int_0^1L\,\big\lVert t(y-x)\big\rVert\,\lVert y-x\rVert\,dt=L\lVert y-x\rVert^2\int_0^1t\,dt
    LL-smoothness at the points x+t(y−x)x+t(y-x) and xx, whose distance is t∥y−x∥t\lVert y-x\rVert for t∈[0,1]t\in[0,1].
  6. Step 6:
    =L2∥y−x∥2=\frac L2\lVert y-x\rVert^2
    ∫01t dt=12\int_0^1t\,dt=\tfrac12. Moving f(x)+∇f(x)⊤(y−x)f(x)+\nabla f(x)^\top(y-x) to the right-hand side gives the claim.

∎

Corollary

Sufficient decrease of one gradient step

If ff is LL-smooth and x+=x−η∇f(x)x^+=x-\eta\nabla f(x), then f(x+)≤f(x)−η(1−Lη2)∥∇f(x)∥2.f(x^+)\le f(x)-\eta\Big(1-\frac{L\eta}2\Big)\lVert\nabla f(x)\rVert^2 . Every step 0<η<2/L0<\eta<2/L decreases ff unless ∇f(x)=0\nabla f(x)=0, and the step η=1/L\eta=1/L maximises the guaranteed decrease: f(x+)≤f(x)−12L∥∇f(x)∥2f(x^+)\le f(x)-\tfrac1{2L}\lVert\nabla f(x)\rVert^2.

Proof (5 steps)
  1. Step 1:
    f(x+)≤f(x)+∇f(x)⊤(−η∇f(x))+L2∥−η∇f(x)∥2f(x^+)\le f(x)+\nabla f(x)^\top\big(-\eta\nabla f(x)\big)+\frac L2\big\lVert-\eta\nabla f(x)\big\rVert^2
    Descent lemma with y=x+y=x^+, so y−x=−η∇f(x)y-x=-\eta\nabla f(x).
  2. Step 2:
    =f(x)−η∥∇f(x)∥2+Lη22∥∇f(x)∥2=f(x)-\eta\lVert\nabla f(x)\rVert^2+\frac{L\eta^2}2\lVert\nabla f(x)\rVert^2
    g⊤g=∥g∥2g^\top g=\lVert g\rVert^2 and ∥−ηg∥2=η2∥g∥2\lVert-\eta g\rVert^2=\eta^2\lVert g\rVert^2.
  3. Step 3:
    =f(x)−η(1−Lη2)∥∇f(x)∥2=f(x)-\eta\Big(1-\frac{L\eta}2\Big)\lVert\nabla f(x)\rVert^2
    Factor out η∥∇f(x)∥2\eta\lVert\nabla f(x)\rVert^2. The bracket is positive exactly when 0<η<2/L0<\eta<2/L.
  4. Step 4:
    ddη(η−Lη22)=1−Lη=0  ⟺  η=1L\frac{d}{d\eta}\Big(\eta-\frac{L\eta^2}2\Big)=1-L\eta=0\iff\eta=\frac1L
    The guaranteed decrease η−L2η2\eta-\tfrac{L}{2}\eta^2 is a concave parabola in η\eta, maximal where its derivative vanishes.
  5. Step 5:
    η=1L:f(x+)≤f(x)−1L(1−12)∥∇f(x)∥2=f(x)−12L∥∇f(x)∥2\eta=\frac1L:\qquad f(x^+)\le f(x)-\frac1L\Big(1-\frac12\Big)\lVert\nabla f(x)\rVert^2=f(x)-\frac1{2L}\lVert\nabla f(x)\rVert^2
    Substitute η=1/L\eta=1/L, so Lη/2=1/2L\eta/2=1/2.

∎

Convergence of gradient descent: the O(1/k) rate and the step-size threshold#

Lemma

First-order characterisation of convexity

If ff is convex and differentiable, i.e. f((1−t)x+ty)≤(1−t)f(x)+tf(y)f\big((1-t)x+ty\big)\le(1-t)f(x)+tf(y) for all x,yx,y and t∈[0,1]t\in[0,1], then f(y)≥f(x)+∇f(x)⊤(y−x)f(y)\ge f(x)+\nabla f(x)^\top(y-x) for all x,yx,y: the tangent plane lies below the graph.

Proof (3 steps)
  1. Step 1:
    f(x+t(y−x))≤f(x)+t(f(y)−f(x)),t∈(0,1]f\big(x+t(y-x)\big)\le f(x)+t\big(f(y)-f(x)\big),\qquad t\in(0,1]
    Rewrite (1−t)x+ty=x+t(y−x)(1-t)x+ty=x+t(y-x) and (1−t)f(x)+tf(y)=f(x)+t(f(y)−f(x))(1-t)f(x)+tf(y)=f(x)+t(f(y)-f(x)) in the definition of convexity.
  2. Step 2:
    f(x+t(y−x))−f(x)t≤f(y)−f(x)\frac{f\big(x+t(y-x)\big)-f(x)}t\le f(y)-f(x)
    Subtract f(x)f(x) and divide by t>0t>0.
  3. Step 3:
    lim⁡t→0+f(x+t(y−x))−f(x)t=∇f(x)⊤(y−x)≤f(y)−f(x)\lim_{t\to0^+}\frac{f\big(x+t(y-x)\big)-f(x)}t=\nabla f(x)^\top(y-x)\le f(y)-f(x)
    The left side is the directional derivative of ff at xx in direction y−xy-x; limits preserve the non-strict inequality.

∎

Theorem

Gradient descent on convex L-smooth functions converges at rate O(1/k)

Let ff be convex and LL-smooth with a minimiser x∗x^\ast, f∗=f(x∗)f^\ast=f(x^\ast), and run xk+1=xk−1L∇f(xk)x_{k+1}=x_k-\tfrac1L\nabla f(x_k) from any x0x_0. Then for every k≥1k\ge1 f(xk)−f∗≤L ∥x0−x∗∥22k,f(x_k)-f^\ast\le\frac{L\,\lVert x_0-x^\ast\rVert^2}{2k}, and the distance ∥xk−x∗∥\lVert x_k-x^\ast\rVert never increases. To reach accuracy ϵ\epsilon it suffices to take k≥L∥x0−x∗∥2/(2ϵ)k\ge L\lVert x_0-x^\ast\rVert^2/(2\epsilon) steps.

Proof (10 steps)
  1. Step 1:
    f(xk)≤f∗+gk⊤(xk−x∗),gk:=∇f(xk)f(x_k)\le f^\ast+g_k^\top(x_k-x^\ast),\qquad g_k:=\nabla f(x_k)
    First-order convexity with x=xkx=x_k, y=x∗y=x^\ast: f∗≥f(xk)+gk⊤(x∗−xk)f^\ast\ge f(x_k)+g_k^\top(x^\ast-x_k); rearrange.
  2. Step 2:
    f(xk+1)≤f(xk)−12L∥gk∥2f(x_{k+1})\le f(x_k)-\frac1{2L}\lVert g_k\rVert^2
    Sufficient decrease with step 1/L1/L (corollary of the descent lemma).
  3. Step 3:
    f(xk+1)−f∗≤gk⊤(xk−x∗)−12L∥gk∥2f(x_{k+1})-f^\ast\le g_k^\top(x_k-x^\ast)-\frac1{2L}\lVert g_k\rVert^2
    Insert the first inequality into the second.
  4. Step 4:
    ∥xk+1−x∗∥2=∥(xk−x∗)−1Lgk∥2=∥xk−x∗∥2−2Lgk⊤(xk−x∗)+1L2∥gk∥2\lVert x_{k+1}-x^\ast\rVert^2=\Big\lVert(x_k-x^\ast)-\tfrac1Lg_k\Big\rVert^2=\lVert x_k-x^\ast\rVert^2-\frac2Lg_k^\top(x_k-x^\ast)+\frac1{L^2}\lVert g_k\rVert^2
    Substitute the update xk+1=xk−gk/Lx_{k+1}=x_k-g_k/L and expand ∥u−v∥2=∥u∥2−2u⊤v+∥v∥2\lVert u-v\rVert^2=\lVert u\rVert^2-2u^\top v+\lVert v\rVert^2.
  5. Step 5:
    gk⊤(xk−x∗)−12L∥gk∥2=L2(∥xk−x∗∥2−∥xk+1−x∗∥2)g_k^\top(x_k-x^\ast)-\frac1{2L}\lVert g_k\rVert^2=\frac L2\Big(\lVert x_k-x^\ast\rVert^2-\lVert x_{k+1}-x^\ast\rVert^2\Big)
    Multiply the expansion by L/2L/2 and solve for gk⊤(xk−x∗)−12L∥gk∥2g_k^\top(x_k-x^\ast)-\tfrac1{2L}\lVert g_k\rVert^2.
  6. Step 6:
    0≤f(xk+1)−f∗≤L2(∥xk−x∗∥2−∥xk+1−x∗∥2)0\le f(x_{k+1})-f^\ast\le\frac L2\Big(\lVert x_k-x^\ast\rVert^2-\lVert x_{k+1}-x^\ast\rVert^2\Big)
    Combine the two previous lines; the left inequality holds because f∗f^\ast is the minimum. In particular ∥xk+1−x∗∥≤∥xk−x∗∥\lVert x_{k+1}-x^\ast\rVert\le\lVert x_k-x^\ast\rVert.
  7. Step 7:
    ∑i=1k(f(xi)−f∗)≤L2∑i=0k−1(∥xi−x∗∥2−∥xi+1−x∗∥2)=L2(∥x0−x∗∥2−∥xk−x∗∥2)\sum_{i=1}^{k}\big(f(x_i)-f^\ast\big)\le\frac L2\sum_{i=0}^{k-1}\Big(\lVert x_i-x^\ast\rVert^2-\lVert x_{i+1}-x^\ast\rVert^2\Big)=\frac L2\Big(\lVert x_0-x^\ast\rVert^2-\lVert x_k-x^\ast\rVert^2\Big)
    Sum the per-step bound over i=0,…,k−1i=0,\dots,k-1; the right-hand side telescopes, leaving the first and the last term.
  8. Step 8:
    ∑i=1k(f(xi)−f∗)≤L2∥x0−x∗∥2\sum_{i=1}^{k}\big(f(x_i)-f^\ast\big)\le\frac L2\lVert x_0-x^\ast\rVert^2
    Drop the non-negative term ∥xk−x∗∥2\lVert x_k-x^\ast\rVert^2.
  9. Step 9:
    f(xk)−f∗≤f(xi)−f∗(1≤i≤k)⟹k(f(xk)−f∗)≤∑i=1k(f(xi)−f∗)f(x_k)-f^\ast\le f(x_i)-f^\ast\quad(1\le i\le k)\quad\Longrightarrow\quad k\big(f(x_k)-f^\ast\big)\le\sum_{i=1}^k\big(f(x_i)-f^\ast\big)
    Monotonicity: by sufficient decrease f(xi+1)≤f(xi)f(x_{i+1})\le f(x_i), so the last iterate has the smallest value; each of the kk summands is at least f(xk)−f∗f(x_k)-f^\ast.
  10. Step 10:
    f(xk)−f∗≤L∥x0−x∗∥22kf(x_k)-f^\ast\le\frac{L\lVert x_0-x^\ast\rVert^2}{2k}
    Chain the two previous inequalities and divide by k≥1k\ge1.

∎

Theorem

Gradient descent on a quadratic: convergence if and only if 0 < η < 2/λ_max

Let f(w)=12(w−w∗)⊤A(w−w∗)f(w)=\tfrac12(w-w^\ast)^\top A(w-w^\ast) with AA symmetric positive definite, eigenvalues 0<λ1≤⋯≤λn=λmax⁡0<\lambda_1\le\dots\le\lambda_n=\lambda_{\max} and orthonormal eigenvectors q1,…,qnq_1,\dots,q_n. Gradient descent with constant step η\eta satisfies f(wk)−f∗=12∑j=1nλj (1−ηλj)2k u0,j2,u0,j=qj⊤(w0−w∗).f(w_k)-f^\ast=\frac12\sum_{j=1}^n\lambda_j\,(1-\eta\lambda_j)^{2k}\,u_{0,j}^2,\qquad u_{0,j}=q_j^\top(w_0-w^\ast). Hence f(wk)→f∗=0f(w_k)\to f^\ast=0 for every start if and only if 0<η<2/λmax⁡0<\eta<2/\lambda_{\max}, in which case f(wk)−f∗≤ρ2k(f(w0)−f∗)f(w_k)-f^\ast\le\rho^{2k}(f(w_0)-f^\ast) with ρ=max⁡j∣1−ηλj∣<1\rho=\max_j\lvert1-\eta\lambda_j\rvert<1; if η>2/λmax⁡\eta>2/\lambda_{\max} and u0,n≠0u_{0,n}\ne0, then f(wk)→∞f(w_k)\to\infty.

Proof (9 steps)
  1. Step 1:
    ∇f(w)=A(w−w∗)\nabla f(w)=A(w-w^\ast)
    Gradient of the quadratic form 12v⊤Av\tfrac12v^\top Av is AvAv for symmetric AA, here with v=w−w∗v=w-w^\ast.
  2. Step 2:
    wk+1−w∗=(wk−w∗)−ηA(wk−w∗)=(I−ηA)(wk−w∗)w_{k+1}-w^\ast=(w_k-w^\ast)-\eta A(w_k-w^\ast)=(I-\eta A)(w_k-w^\ast)
    Subtract w∗w^\ast from both sides of the GD update and factor.
  3. Step 3:
    wk−w∗=(I−ηA)k(w0−w∗)w_k-w^\ast=(I-\eta A)^k(w_0-w^\ast)
    Induction on kk with the one-step recursion.
  4. Step 4:
    A=QΛQ⊤⟹(I−ηA)k=Q (I−ηΛ)k Q⊤A=Q\Lambda Q^\top\quad\Longrightarrow\quad(I-\eta A)^k=Q\,(I-\eta\Lambda)^k\,Q^\top
    Spectral theorem for the real symmetric AA (QQ orthogonal with columns qjq_j); I−ηA=Q(I−ηΛ)Q⊤I-\eta A=Q(I-\eta\Lambda)Q^\top and Q⊤Q=IQ^\top Q=I makes the inner factors cancel in the power.
  5. Step 5:
    uk:=Q⊤(wk−w∗)⟹uk,j=(1−ηλj)k u0,ju_k:=Q^\top(w_k-w^\ast)\quad\Longrightarrow\quad u_{k,j}=(1-\eta\lambda_j)^k\,u_{0,j}
    Multiply the previous identity by Q⊤Q^\top; (I−ηΛ)k(I-\eta\Lambda)^k is diagonal, so the coordinates along the eigenvectors decouple.
  6. Step 6:
    f(wk)−f∗=12uk⊤Λuk=12∑j=1nλj(1−ηλj)2ku0,j2f(w_k)-f^\ast=\frac12u_k^\top\Lambda u_k=\frac12\sum_{j=1}^n\lambda_j(1-\eta\lambda_j)^{2k}u_{0,j}^2
    (wk−w∗)⊤A(wk−w∗)=uk⊤Q⊤AQuk=uk⊤Λuk(w_k-w^\ast)^\top A(w_k-w^\ast)=u_k^\top Q^\top AQu_k=u_k^\top\Lambda u_k, and f∗=f(w∗)=0f^\ast=f(w^\ast)=0.
  7. Step 7:
    0<η<2λmax⁡  ⟺  −1<1−ηλj<1  ∀j⟹f(wk)−f∗≤ρ2k(f(w0)−f∗)→00<\eta<\frac2{\lambda_{\max}}\iff-1<1-\eta\lambda_j<1\ \ \forall j\quad\Longrightarrow\quad f(w_k)-f^\ast\le\rho^{2k}\big(f(w_0)-f^\ast\big)\to0
    For λj∈(0,λmax⁡]\lambda_j\in(0,\lambda_{\max}]: 1−ηλj<11-\eta\lambda_j<1 iff η>0\eta>0, and 1−ηλj>−11-\eta\lambda_j>-1 for all jj iff ηλmax⁡<2\eta\lambda_{\max}<2. Then each factor (1−ηλj)2k≤ρ2k(1-\eta\lambda_j)^{2k}\le\rho^{2k} with ρ<1\rho<1.
  8. Step 8:
    η≥2λmax⁡ or η≤0:∣1−ηλn∣≥1⟹f(wk)−f∗≥12λnu0,n2>0  ∀k\eta\ge\frac2{\lambda_{\max}}\ \text{or}\ \eta\le0:\quad\lvert1-\eta\lambda_n\rvert\ge1\quad\Longrightarrow\quad f(w_k)-f^\ast\ge\frac12\lambda_n u_{0,n}^2>0\ \ \forall k
    Keep only the j=nj=n term of the non-negative sum. If u0,n≠0u_{0,n}\ne0 the error never drops below a positive constant, so convergence fails for such starts.
  9. Step 9:
    η>2λmax⁡:∣1−ηλn∣=ηλn−1>1⟹f(wk)−f∗≥12λn(ηλn−1)2ku0,n2→∞\eta>\frac2{\lambda_{\max}}:\quad\lvert1-\eta\lambda_n\rvert=\eta\lambda_n-1>1\quad\Longrightarrow\quad f(w_k)-f^\ast\ge\frac12\lambda_n(\eta\lambda_n-1)^{2k}u_{0,n}^2\to\infty
    1−ηλn<−11-\eta\lambda_n<-1, so the qnq_n-coordinate flips sign and grows geometrically: the iterates oscillate with increasing amplitude.

∎

Derivation

Example: f(w) = ½(w₁² + 10w₂²) from w₀ = (−1.8, 1.2)

  1. Step 1:
    A=(10010),λmax⁡=10,2λmax⁡=0.2,wk=(−1.8 (1−η)k, 1.2 (1−10η)k)A=\begin{pmatrix}1&0\\0&10\end{pmatrix},\quad\lambda_{\max}=10,\quad\frac2{\lambda_{\max}}=0.2,\quad w_k=\big(-1.8\,(1-\eta)^k,\ 1.2\,(1-10\eta)^k\big)
    AA is already diagonal, so Q=IQ=I and each coordinate evolves with its own factor 1−ηλj1-\eta\lambda_j; w∗=0w^\ast=0.
  2. Step 2:
    η=0.15:1−η=0.85,1−10η=−0.5,ρ=0.85\eta=0.15:\quad1-\eta=0.85,\quad1-10\eta=-0.5,\quad\rho=0.85
    Both factors lie in (−1,1)(-1,1): w2w_2 alternates in sign and halves every step, w1w_1 shrinks slowly. After 4040 steps ∣w1∣=1.8⋅0.8540≈2.7×10−3\lvert w_1\rvert=1.8\cdot0.85^{40}\approx2.7\times10^{-3}.
  3. Step 3:
    η=0.21:1−10η=−1.1,∣wk,2∣=1.2⋅1.1k\eta=0.21:\quad1-10\eta=-1.1,\quad\lvert w_{k,2}\rvert=1.2\cdot1.1^k
    0.21>0.2=2/λmax⁡0.21>0.2=2/\lambda_{\max}: the steep coordinate overshoots further at every step. After 2525 steps ∣w2∣=1.2⋅1.125≈13.0\lvert w_2\rvert=1.2\cdot1.1^{25}\approx13.0, although the w1w_1 factor 0.790.79 still contracts.
  4. Step 4:
    ρ(η)=max⁡(∣1−η∣,∣1−10η∣) is minimal at 1−η=10η−1  ⟺  η=211,ρ=911\rho(\eta)=\max\big(\lvert1-\eta\rvert,\lvert1-10\eta\rvert\big)\ \text{is minimal at}\ 1-\eta=10\eta-1\iff\eta=\frac2{11},\quad\rho=\frac9{11}
    On (0,0.2)(0,0.2), ∣1−η∣=1−η\lvert1-\eta\rvert=1-\eta decreases while ∣1−10η∣\lvert1-10\eta\rvert decreases until η=0.1\eta=0.1 and then increases; the maximum of the two is smallest where the decreasing and the increasing branch cross. In general η=2/(λmin⁡+λmax⁡)\eta=2/(\lambda_{\min}+\lambda_{\max}) gives ρ=(κ−1)/(κ+1)\rho=(\kappa-1)/(\kappa+1) with κ=λmax⁡/λmin⁡\kappa=\lambda_{\max}/\lambda_{\min}.

The first plot uses η=0.15<0.2\eta=0.15<0.2: the trajectory zig-zags across the narrow valley along w2w_2 while creeping along the flat direction w1w_1. The second plot raises the step to η=0.21>2/λmax⁡\eta=0.21>2/\lambda_{\max} with everything else unchanged; drag the learning-rate slider across 0.20.2 to watch the transition from damped to growing oscillations.

Quick check

Gradient descent minimises f(w)=12(3w12+25w22)f(w)=\tfrac12\big(3w_1^2+25w_2^2\big). What is the supremum of the constant step sizes η\eta for which it converges from every starting point?

tolerance ±1e-9

Stochastic gradients and the bridge to shot noise#

For a training loss L(θ)=1m∑i=1mℓi(θ)\mathcal L(\theta)=\tfrac1m\sum_{i=1}^m\ell_i(\theta), one exact gradient costs mm per-sample gradients. Stochastic gradient descent (SGD) replaces it by a cheap random estimate: draw indices I1,…,IBI_1,\dots,I_B independently and uniformly from {1,…,m}\{1,\dots,m\} (sampling with replacement) and use the minibatch gradient g^(θ)=1B∑b=1B∇ℓIb(θ),θk+1=θk−η g^(θk),\hat g(\theta)=\frac1B\sum_{b=1}^B\nabla\ell_{I_b}(\theta),\qquad\theta_{k+1}=\theta_k-\eta\,\hat g(\theta_k), with fresh indices at every step. Write Σ(θ)=1m∑i=1m(∇ℓi−∇L)(∇ℓi−∇L)⊤\Sigma(\theta)=\tfrac1m\sum_{i=1}^m\big(\nabla\ell_i-\nabla\mathcal L\big)\big(\nabla\ell_i-\nabla\mathcal L\big)^\top for the covariance of one uniformly drawn per-sample gradient.

Theorem

The minibatch gradient is unbiased with covariance Σ/B

For fixed θ\theta, E[g^(θ)]=∇L(θ),Cov(g^(θ))=Σ(θ)B,E∥g^(θ)−∇L(θ)∥2=tr⁡Σ(θ)B.\mathbb E\big[\hat g(\theta)\big]=\nabla\mathcal L(\theta),\qquad\mathrm{Cov}\big(\hat g(\theta)\big)=\frac{\Sigma(\theta)}B,\qquad\mathbb E\big\lVert\hat g(\theta)-\nabla\mathcal L(\theta)\big\rVert^2=\frac{\operatorname{tr}\Sigma(\theta)}B .

Proof (6 steps)
  1. Step 1:
    E[∇ℓIb]=∑i=1mPr⁡(Ib=i) ∇ℓi=1m∑i=1m∇ℓi=∇L\mathbb E\big[\nabla\ell_{I_b}\big]=\sum_{i=1}^m\Pr(I_b=i)\,\nabla\ell_i=\frac1m\sum_{i=1}^m\nabla\ell_i=\nabla\mathcal L
    Expectation of a function of the discrete uniform index IbI_b, with Pr⁡(Ib=i)=1m\Pr(I_b=i)=\tfrac1m; the last equality is linearity of the gradient.
  2. Step 2:
    E[g^]=1B∑b=1BE[∇ℓIb]=1B⋅B ∇L=∇L\mathbb E[\hat g]=\frac1B\sum_{b=1}^B\mathbb E\big[\nabla\ell_{I_b}\big]=\frac1B\cdot B\,\nabla\mathcal L=\nabla\mathcal L
    Linearity of expectation; every index has the same distribution. This is unbiasedness.
  3. Step 3:
    g^−∇L=1B∑b=1Bξb,ξb:=∇ℓIb−∇L,E[ξb]=0\hat g-\nabla\mathcal L=\frac1B\sum_{b=1}^B\xi_b,\qquad\xi_b:=\nabla\ell_{I_b}-\nabla\mathcal L,\qquad\mathbb E[\xi_b]=0
    Subtract ∇L=1B∑b∇L\nabla\mathcal L=\tfrac1B\sum_b\nabla\mathcal L; the ξb\xi_b are independent because the indices are, and centred by the first step.
  4. Step 4:
    Cov(g^)=E[1B2∑b,b′ξbξb′⊤]=1B2∑b=1BE[ξbξb⊤]\mathrm{Cov}(\hat g)=\mathbb E\Big[\frac1{B^2}\sum_{b,b'}\xi_b\xi_{b'}^\top\Big]=\frac1{B^2}\sum_{b=1}^B\mathbb E\big[\xi_b\xi_b^\top\big]
    For b≠b′b\ne b', independence and zero mean give E[ξbξb′⊤]=E[ξb]E[ξb′]⊤=0\mathbb E[\xi_b\xi_{b'}^\top]=\mathbb E[\xi_b]\mathbb E[\xi_{b'}]^\top=0; only the BB diagonal terms remain.
  5. Step 5:
    E[ξbξb⊤]=1m∑i=1m(∇ℓi−∇L)(∇ℓi−∇L)⊤=Σ\mathbb E\big[\xi_b\xi_b^\top\big]=\frac1m\sum_{i=1}^m\big(\nabla\ell_i-\nabla\mathcal L\big)\big(\nabla\ell_i-\nabla\mathcal L\big)^\top=\Sigma
    Same uniform average as in the first step, applied to the outer product.
  6. Step 6:
    Cov(g^)=B ΣB2=ΣB,E∥g^−∇L∥2=tr⁡Cov(g^)=tr⁡ΣB\mathrm{Cov}(\hat g)=\frac{B\,\Sigma}{B^2}=\frac\Sigma B,\qquad\mathbb E\lVert\hat g-\nabla\mathcal L\rVert^2=\operatorname{tr}\mathrm{Cov}(\hat g)=\frac{\operatorname{tr}\Sigma}B
    Sum the BB identical terms; the mean squared error of a vector estimator is the trace of its covariance, since ∥v∥2=tr⁡(vv⊤)\lVert v\rVert^2=\operatorname{tr}(vv^\top).

∎

Proposition

Expected progress of one SGD step

If L\mathcal L is LL-smooth, then one SGD step from a fixed θ\theta satisfies E[L(θ−ηg^)]≤L(θ)−η(1−Lη2)∥∇L(θ)∥2+Lη22 tr⁡Σ(θ)B.\mathbb E\big[\mathcal L(\theta-\eta\hat g)\big]\le\mathcal L(\theta)-\eta\Big(1-\frac{L\eta}2\Big)\lVert\nabla\mathcal L(\theta)\rVert^2+\frac{L\eta^2}{2}\,\frac{\operatorname{tr}\Sigma(\theta)}B . The last term does not vanish at a minimiser: with a constant step, SGD stalls in a noise ball whose size shrinks with smaller η\eta or larger BB.

Proof (5 steps)
  1. Step 1:
    L(θ−ηg^)≤L(θ)−η ∇L(θ)⊤g^+Lη22∥g^∥2\mathcal L(\theta-\eta\hat g)\le\mathcal L(\theta)-\eta\,\nabla\mathcal L(\theta)^\top\hat g+\frac{L\eta^2}2\lVert\hat g\rVert^2
    Descent lemma with y=θ−ηg^y=\theta-\eta\hat g, valid for every realisation of the minibatch.
  2. Step 2:
    E[∇L(θ)⊤g^]=∇L(θ)⊤E[g^]=∥∇L(θ)∥2\mathbb E\big[\nabla\mathcal L(\theta)^\top\hat g\big]=\nabla\mathcal L(\theta)^\top\mathbb E[\hat g]=\lVert\nabla\mathcal L(\theta)\rVert^2
    θ\theta is fixed, so ∇L(θ)\nabla\mathcal L(\theta) is a constant vector; then use unbiasedness.
  3. Step 3:
    E∥g^∥2=∥Eg^∥2+E∥g^−Eg^∥2=∥∇L∥2+tr⁡ΣB\mathbb E\lVert\hat g\rVert^2=\big\lVert\mathbb E\hat g\big\rVert^2+\mathbb E\big\lVert\hat g-\mathbb E\hat g\big\rVert^2=\lVert\nabla\mathcal L\rVert^2+\frac{\operatorname{tr}\Sigma}B
    Write g^=Eg^+(g^−Eg^)\hat g=\mathbb E\hat g+(\hat g-\mathbb E\hat g) and expand the square; the cross term has expectation 00. Then apply the theorem above.
  4. Step 4:
    E[L(θ−ηg^)]≤L(θ)−η∥∇L∥2+Lη22(∥∇L∥2+tr⁡ΣB)\mathbb E\big[\mathcal L(\theta-\eta\hat g)\big]\le\mathcal L(\theta)-\eta\lVert\nabla\mathcal L\rVert^2+\frac{L\eta^2}2\Big(\lVert\nabla\mathcal L\rVert^2+\frac{\operatorname{tr}\Sigma}B\Big)
    Take expectations of the first inequality and substitute the two previous lines.
  5. Step 5:
    =L(θ)−η(1−Lη2)∥∇L(θ)∥2+Lη22 tr⁡Σ(θ)B=\mathcal L(\theta)-\eta\Big(1-\frac{L\eta}2\Big)\lVert\nabla\mathcal L(\theta)\rVert^2+\frac{L\eta^2}{2}\,\frac{\operatorname{tr}\Sigma(\theta)}B
    Collect the two terms in ∥∇L∥2\lVert\nabla\mathcal L\rVert^2.

∎

The plot runs the same quadratic with step 0.10.1 and zero-mean noise added to every gradient: the trajectory reaches the valley as fast as GD but then jitters around the minimiser instead of converging, as the noise term of the proposition predicts. Lowering the noise level or the learning rate shrinks the jitter.

The bridge to quantum gradients. A variational quantum circuit prepares ∣ψ(θ)⟩\lvert\psi(\theta)\rangle and its model output is an expectation value E(θ)=⟨ψ(θ)∣P∣ψ(θ)⟩E(\theta)=\langle\psi(\theta)\rvert P\lvert\psi(\theta)\rangle of an observable, here a Pauli string PP with eigenvalues ±1\pm1. Hardware never returns E(θ)E(\theta) itself, only measurement outcomes, so every quantity entering gradient descent is an estimate — exactly the situation of SGD.

Derivation

Shot noise: unbiased estimates with variance proportional to 1/N

  1. Step 1:
    P=Π+−Π−,Π±=I±P2,Pr⁡(s=±1)=⟨ψ∣Π±∣ψ⟩=1±μ2,μ:=⟨ψ∣P∣ψ⟩P=\Pi_+-\Pi_-,\quad\Pi_\pm=\frac{I\pm P}2,\quad\Pr(s=\pm1)=\langle\psi\rvert\Pi_\pm\lvert\psi\rangle=\frac{1\pm\mu}2,\quad\mu:=\langle\psi\rvert P\lvert\psi\rangle
    P2=IP^2=I, so Π±\Pi_\pm are the projectors onto its ±1\pm1 eigenspaces; the Born rule gives the outcome probabilities of one shot s∈{+1,−1}s\in\{+1,-1\}.
  2. Step 2:
    E[s]=(+1)1+μ2+(−1)1−μ2=μ\mathbb E[s]=(+1)\frac{1+\mu}2+(-1)\frac{1-\mu}2=\mu
    Definition of the expectation of a two-valued random variable.
  3. Step 3:
    Var(s)=E[s2]−E[s]2=1−μ2\mathrm{Var}(s)=\mathbb E[s^2]-\mathbb E[s]^2=1-\mu^2
    s2=1s^2=1 for both outcomes.
  4. Step 4:
    E^N=1N∑n=1Nsn:E[E^N]=μ,Var(E^N)=1−μ2N\hat E_N=\frac1N\sum_{n=1}^Ns_n:\qquad\mathbb E[\hat E_N]=\mu,\qquad\mathrm{Var}(\hat E_N)=\frac{1-\mu^2}N
    The NN shots are independent and identically distributed: linearity gives the mean, and variances of independent terms add, each scaled by 1/N21/N^2.
  5. Step 5:
    g^j=12(E^N(θ+π2ej)−E^N(θ−π2ej)),E[g^j]=12(E(θ+π2ej)−E(θ−π2ej))=∂E∂θj\hat g_j=\frac12\Big(\hat E_N\big(\theta+\tfrac\pi2e_j\big)-\hat E_N\big(\theta-\tfrac\pi2e_j\big)\Big),\qquad\mathbb E[\hat g_j]=\frac12\big(E(\theta+\tfrac\pi2e_j)-E(\theta-\tfrac\pi2e_j)\big)=\frac{\partial E}{\partial\theta_j}
    Linearity of expectation; the last equality is the parameter-shift rule, exact for gates e−iθjG/2e^{-i\theta_jG/2} with a Pauli generator GG, which Stage 4 proves. The shot-based gradient is therefore unbiased.
  6. Step 6:
    Var(g^j)=14(1−μ+2N+1−μ−2N)≤12N\mathrm{Var}(\hat g_j)=\frac14\Big(\frac{1-\mu_+^2}N+\frac{1-\mu_-^2}N\Big)\le\frac1{2N}
    The two shifted circuits are measured with independent shots, so the variances add with the factor (12)2(\tfrac12)^2; each 1−μ±2≤11-\mu_\pm^2\le1.

Result

E[g^j]=∂jE(θ),Var(g^j)≤12N\boxed{\mathbb E[\hat g_j]=\partial_jE(\theta),\qquad\mathrm{Var}(\hat g_j)\le\frac1{2N}}

Quick check

The minibatch size is raised from B=8B=8 to B=32B=32 (sampling with replacement). What happens to the standard deviation of each coordinate of the gradient estimate g^\hat g?

Select one answer

References

  1. T. Hastie, R. Tibshirani and J. Friedman, The Elements of Statistical Learning: Data Mining, Inference, and Prediction, 2nd ed., Springer (2009), §2.9, §3.2, §4.4, §7.3.doi:10.1007/978-0-387-84858-7
  2. Y. Nesterov, Introductory Lectures on Convex Optimization: A Basic Course, Kluwer Academic Publishers (2004), §1.2.2 and §2.1.doi:10.1007/978-1-4419-8853-9
  3. L. Bottou, F. E. Curtis and J. Nocedal, Optimization Methods for Large-Scale Machine Learning, SIAM Review 60(2), 223–311 (2018).doi:10.1137/16M1080173
  4. S. Shalev-Shwartz and S. Ben-David, Understanding Machine Learning: From Theory to Algorithms, Cambridge University Press (2014), ch. 2 and ch. 14.doi:10.1017/CBO9781107298019
  5. M. Schuld, V. Bergholm, C. Gogolin, J. Izaac and N. Killoran, Evaluating analytic gradients on quantum hardware, Physical Review A 99, 032331 (2019).doi:10.1103/PhysRevA.99.032331

Hands-on labs

  • Gradient of the mean-squared error

    1 step · auto-graded

    Implement the vectorised gradient of the mean-squared-error loss L(w,b)=1m∑i=1m(w⊤xi+b−yi)2\mathcal L(w,b)=\tfrac1m\sum_{i=1}^m(w^\top x_i+b-y_i)^2 of a linear model, namely ∇wL=2m∑irixi=2mX⊤r\nabla_w\mathcal L=\tfrac2m\sum_i r_ix_i=\tfrac2mX^\top r and ∂bL=2m∑iri\partial_b\mathcal L=\tfrac2m\sum_ir_i with residuals ri=w⊤xi+b−yir_i=w^\top x_i+b-y_i. The function is graded on one visible worked example and five hidden random problems with up to 5050 samples and 55 features, to absolute tolerance 10−910^{-9}.

    Open lab