QuantumCore

Lesson 48 · 110 min · advanced

The Parameter-Shift Rule

Loading progress…

The parameter-shift rule computes the exact derivative of an expectation value f(θ)=⟨ψ∣U†(θ)BU(θ)∣ψ⟩f(\theta)=\langle\psi\vert U^\dagger(\theta)BU(\theta)\vert\psi\rangle from two further evaluations of the same circuit at shifted parameter values: ∂jf=12[f(θ+π2ej)−f(θ−π2ej)]\partial_jf=\tfrac12\big[f(\theta+\tfrac\pi2e_j)-f(\theta-\tfrac\pi2e_j)\big] for a Pauli rotation. This lesson derives the rule from the spectral decomposition of the gate generator. A generator with eigenvalues ±r\pm r makes ff a single sinusoid of frequency 2r2r, and a sum-to-product identity then gives the general rule f′(θ)=rsin⁡2rs[f(θ+s)−f(θ−s)]f'(\theta)=\frac{r}{\sin2rs}\big[f(\theta+s)-f(\theta-s)\big]. The two-term rule is proved both in the Schrödinger and in the Heisenberg picture. The lesson then shows that the estimator is exact and unbiased under shot noise, extends the rule to parameters that occur in several gates, derives the four-term rule for controlled rotations by solving a 2×22\times2 linear system, and counts the circuit evaluations a gradient costs.

Prerequisites: Why Backpropagation Fails on Quantum Hardware, Spectral Theorem & Functions of Operators

Objectives and key results

  • Reduce the derivative with respect to one gate angle to a function f(θ)=⟨ψ∣U†(θ)BU(θ)∣ψ⟩f(\theta)=\langle\psi\vert U^\dagger(\theta)BU(\theta)\vert\psi\rangle by absorbing the gates before and after the parameterised gate into ∣ψ⟩\lvert\psi\rangle and BB.

  • Prove that a generator with eigenvalues ±r\pm r makes ff a single sinusoid a0+a1cos⁡2rθ+b1sin⁡2rθa_0+a_1\cos2r\theta+b_1\sin2r\theta, with a0,a1,b1a_0,a_1,b_1 expressed through the spectral projectors Π±\Pi_\pm.

  • Derive the general shift rule f′(θ)=rsin⁡2rs[f(θ+s)−f(θ−s)]f'(\theta)=\frac{r}{\sin2rs}\big[f(\theta+s)-f(\theta-s)\big], explain the choice s=π4rs=\frac{\pi}{4r}, and specialise it to the two-term rule for Pauli rotations in both the Schrödinger and the Heisenberg picture.

  • Show that the parameter-shift estimator is exact and unbiased under shot noise, compute its variance, and compare it with the bias and variance of finite differences.

  • Differentiate parameters that occur in several gates with the chain rule, and derive the four-term rule for controlled rotations from the frequencies {12,1}\{\tfrac12,1\} of their generators.

  • Implement a parameter-shift gradient that uses exactly 2P2P circuit evaluations, and count the cost of gradients with controlled rotations.

theorem · Two-term parameter-shift rule (Pauli rotations)
∂⟨H⟩∂θj=12[⟨H⟩θj+π2−⟨H⟩θj−π2],RP(θj)=e−iθjP/2\frac{\partial\langle H\rangle}{\partial\theta_j}=\frac12\Big[\langle H\rangle_{\theta_j+\frac\pi2}-\langle H\rangle_{\theta_j-\frac\pi2}\Big],\qquad R_P(\theta_j)=e^{-i\theta_jP/2}
identity · Single-frequency form
G=r(Π+−Π−) ⇒ f(θ)=a0+a1cos⁡2rθ+b1sin⁡2rθ,a0=⟨Π+BΠ+⟩+⟨Π−BΠ−⟩,a1−ib1=2⟨Π+BΠ−⟩G=r(\Pi_+-\Pi_-)\ \Rightarrow\ f(\theta)=a_0+a_1\cos2r\theta+b_1\sin2r\theta,\quad a_0=\langle\Pi_+B\Pi_+\rangle+\langle\Pi_-B\Pi_-\rangle,\quad a_1-ib_1=2\langle\Pi_+B\Pi_-\rangle
theorem · General shift rule
f(θ+s)−f(θ−s)=sin⁡2rsr f′(θ) ⇒ f′(θ)=rsin⁡2rs[f(θ+s)−f(θ−s)],s⋆=π4rf(\theta+s)-f(\theta-s)=\frac{\sin2rs}{r}\,f'(\theta)\ \Rightarrow\ f'(\theta)=\frac{r}{\sin2rs}\big[f(\theta+s)-f(\theta-s)\big],\qquad s^\star=\frac{\pi}{4r}
identity · Heisenberg form of the rule
ddθU†BU=iU†[G,B]U,i[G,B]=12[U†(π2)BU(π2)−U†(−π2)BU(−π2)](G=P2)\frac{d}{d\theta}U^\dagger BU=iU^\dagger[G,B]U,\qquad i[G,B]=\tfrac12\Big[U^\dagger\big(\tfrac\pi2\big)BU\big(\tfrac\pi2\big)-U^\dagger\big(-\tfrac\pi2\big)BU\big(-\tfrac\pi2\big)\Big]\quad\big(G=\tfrac P2\big)
formula · Four-term rule for controlled rotations
f′(θ)=c+[f(θ+π2)−f(θ−π2)]−c−[f(θ+3π2)−f(θ−3π2)],c±=2±142f'(\theta)=c_+\Big[f\big(\theta+\tfrac\pi2\big)-f\big(\theta-\tfrac\pi2\big)\Big]-c_-\Big[f\big(\theta+\tfrac{3\pi}2\big)-f\big(\theta-\tfrac{3\pi}2\big)\Big],\qquad c_\pm=\frac{\sqrt2\pm1}{4\sqrt2}
formula · Chain rule and cost
dFdα=∑gag ∂f∂θg  (θg=agα+bg),#evaluations=2 mPauli+4 mcontrolled\frac{dF}{d\alpha}=\sum_g a_g\,\frac{\partial f}{\partial\theta_g}\ \ (\theta_g=a_g\alpha+b_g),\qquad\#\text{evaluations}=2\,m_{\text{Pauli}}+4\,m_{\text{controlled}}

The expectation value as a function of one angle#

The previous lesson established two facts. Reverse-mode differentiation is not available on quantum hardware: the intermediate states cannot be read out without tomography, and they cannot be copied. Finite differences are available but poor: the step hh introduces a bias of order h2h^2, and shot noise is amplified by a factor 1/h21/h^2 in the variance. This lesson shows that the derivative of an expectation value with respect to a gate angle can nevertheless be obtained exactly from circuit evaluations, provided the evaluations are taken at large, fixed shifts of the angle instead of infinitesimal ones. The reason is structural: as a function of one angle, an expectation value is a trigonometric polynomial with very few frequencies, and the derivative of such a function is a finite linear combination of its values.

Definition

Gate-isolated cost function

Let W(θ)W(\theta) be a parameterised circuit on nn qubits, HH a Hermitian observable and F(θ)=⟨0n∣W†(θ)HW(θ)∣0n⟩F(\theta)=\langle0^n\vert W^\dagger(\theta)HW(\theta)\vert0^n\rangle. Fix all parameters except the angle θj\theta_j of one gate U(θj)=e−iθjGU(\theta_j)=e^{-i\theta_jG} with Hermitian generator GG, and write W=V>U(θj)V<W=V_>U(\theta_j)V_<, where V<V_< collects the gates before and V>V_> the gates after it. The gate-isolated cost function is f(θ)=⟨ψ∣U†(θ) B U(θ)∣ψ⟩,∣ψ⟩=V<∣0n⟩,B=V>†HV>.f(\theta)=\langle\psi\vert U^\dagger(\theta)\,B\,U(\theta)\vert\psi\rangle,\qquad\lvert\psi\rangle=V_<\lvert0^n\rangle,\qquad B=V_>^\dagger HV_> . For a Pauli rotation RP(θ)=e−iθP/2R_P(\theta)=e^{-i\theta P/2} the generator is G=P/2G=P/2.

Derivation

Absorbing the rest of the circuit into the state and the observable

  1. Step 1:
    W(θ)=V> U(θj) V<,∣Ψ(θ)⟩=W(θ)∣0n⟩W(\theta)=V_>\,U(\theta_j)\,V_<,\qquad\lvert\Psi(\theta)\rangle=W(\theta)\lvert0^n\rangle
    A circuit is a product of gates in time order; the gate carrying θj\theta_j splits it into a part before and a part after. Neither V<V_< nor V>V_> depends on θj\theta_j when θj\theta_j occurs in one gate only.
  2. Step 2:
    F(θ)=⟨0n∣V<† U†(θj) V>† H V> U(θj) V<∣0n⟩F(\theta)=\langle0^n\vert V_<^\dagger\,U^\dagger(\theta_j)\,V_>^\dagger\,H\,V_>\,U(\theta_j)\,V_<\vert0^n\rangle
    The adjoint of a product is the product of the adjoints in reverse order: (V>UV<)†=V<†U†V>†(V_>UV_<)^\dagger=V_<^\dagger U^\dagger V_>^\dagger.
  3. Step 3:
    F(θ)=⟨ψ∣U†(θj) B U(θj)∣ψ⟩=f(θj),∣ψ⟩=V<∣0n⟩,  B=V>†HV>F(\theta)=\langle\psi\vert U^\dagger(\theta_j)\,B\,U(\theta_j)\vert\psi\rangle=f(\theta_j),\qquad\lvert\psi\rangle=V_<\lvert0^n\rangle,\ \ B=V_>^\dagger HV_>
    Group the factors: the gates before the parameterised gate act on the initial state, the gates after it act on the observable.
  4. Step 4:
    B†=V>†H†V>=B,spec⁡(B)=spec⁡(H)B^\dagger=V_>^\dagger H^\dagger V_>=B,\qquad\operatorname{spec}(B)=\operatorname{spec}(H)
    HH is Hermitian, and conjugation by the unitary V>V_> preserves Hermiticity and eigenvalues. In particular ∥B∥=∥H∥\lVert B\rVert=\lVert H\rVert, so ∣f∣≤∥H∥\lvert f\rvert\le\lVert H\rVert.
  5. Step 5:
    U(θ)=e−iθG,G=G†;RP(θ)=e−iθP/2 ⇒ G=12PU(\theta)=e^{-i\theta G},\quad G=G^\dagger;\qquad R_P(\theta)=e^{-i\theta P/2}\ \Rightarrow\ G=\tfrac12P
    Every one-parameter gate of the platform has this form. For a Pauli string PP the generator P/2P/2 is Hermitian with P2=IP^2=I.
  6. Step 6:
    ∂F∂θj(θ)=f′(θj)\frac{\partial F}{\partial\theta_j}(\theta)=f'(\theta_j)
    All other parameters are frozen inside ∣ψ⟩\lvert\psi\rangle and BB, so the partial derivative of the full cost function is the ordinary derivative of the one-variable function ff.

Result

f(θ)=⟨ψ∣U†(θ) B U(θ)∣ψ⟩,U(θ)=e−iθG\boxed{f(\theta)=\langle\psi\vert U^\dagger(\theta)\,B\,U(\theta)\vert\psi\rangle,\qquad U(\theta)=e^{-i\theta G}}

Generators with two eigenvalues#

The dependence of ff on θ\theta is fixed by the spectrum of the generator. The spectral theorem of Stage 1 writes a Hermitian operator as a weighted sum of orthogonal projectors, and an operator function such as e−iθGe^{-i\theta G} acts on each projector separately. When GG has only two eigenvalues, this decomposition makes ff a single sinusoid.

Theorem

Expectation values of two-eigenvalue gates are single sinusoids

Let GG be Hermitian with exactly the two eigenvalues ±r\pm r, r>0r>0, and let Π±\Pi_\pm be the orthogonal projectors onto the eigenspaces, so G=r(Π+−Π−)G=r(\Pi_+-\Pi_-). Then for every state ∣ψ⟩\lvert\psi\rangle and Hermitian BB, f(θ)=a0+a1cos⁡2rθ+b1sin⁡2rθ,f(\theta)=a_0+a_1\cos2r\theta+b_1\sin2r\theta, a0=⟨ψ∣Π+BΠ+∣ψ⟩+⟨ψ∣Π−BΠ−∣ψ⟩,a1=2 Re ⟨ψ∣Π+BΠ−∣ψ⟩,b1=−2 Im ⟨ψ∣Π+BΠ−∣ψ⟩,a_0=\langle\psi\vert\Pi_+B\Pi_+\vert\psi\rangle+\langle\psi\vert\Pi_-B\Pi_-\vert\psi\rangle,\qquad a_1=2\,\mathrm{Re}\,\langle\psi\vert\Pi_+B\Pi_-\vert\psi\rangle,\qquad b_1=-2\,\mathrm{Im}\,\langle\psi\vert\Pi_+B\Pi_-\vert\psi\rangle, all three real. For a Pauli rotation, r=12r=\tfrac12, Π±=12(I±P)\Pi_\pm=\tfrac12(I\pm P) and f(θ)=a0+a1cos⁡θ+b1sin⁡θf(\theta)=a_0+a_1\cos\theta+b_1\sin\theta has period 2π2\pi.

Proof (8 steps)
  1. Step 1:
    G=rΠ+−rΠ−,Π++Π−=I,Π±2=Π±=Π±†,Π+Π−=Π−Π+=0G=r\Pi_+-r\Pi_-,\qquad\Pi_++\Pi_-=I,\qquad\Pi_\pm^2=\Pi_\pm=\Pi_\pm^\dagger,\qquad\Pi_+\Pi_-=\Pi_-\Pi_+=0
    Spectral theorem for the Hermitian operator GG: the eigenspaces of distinct eigenvalues are orthogonal and together span the space.
  2. Step 2:
    e−iθG=e−irθ Π++eirθ Π−,U†(θ)=eirθ Π++e−irθ Π−e^{-i\theta G}=e^{-ir\theta}\,\Pi_++e^{ir\theta}\,\Pi_-,\qquad U^\dagger(\theta)=e^{ir\theta}\,\Pi_++e^{-ir\theta}\,\Pi_-
    A function of a Hermitian operator applies the function to each eigenvalue: h(G)=∑λh(λ)Πλh(G)=\sum_\lambda h(\lambda)\Pi_\lambda, here with h(λ)=e−iθλh(\lambda)=e^{-i\theta\lambda}. The adjoint conjugates the scalar coefficients.
  3. Step 3:
    U†BU=Π+BΠ++Π−BΠ−+e2irθ Π+BΠ−+e−2irθ Π−BΠ+U^\dagger BU=\Pi_+B\Pi_++\Pi_-B\Pi_-+e^{2ir\theta}\,\Pi_+B\Pi_-+e^{-2ir\theta}\,\Pi_-B\Pi_+
    Multiply out the four products. The phases combine as eirθe−irθ=1e^{ir\theta}e^{-ir\theta}=1 on the diagonal blocks and as e±irθe±irθ=e±2irθe^{\pm ir\theta}e^{\pm ir\theta}=e^{\pm2ir\theta} on the off-diagonal blocks.
  4. Step 4:
    z:=⟨ψ∣Π+BΠ−∣ψ⟩,⟨ψ∣Π−BΠ+∣ψ⟩=z‾z:=\langle\psi\vert\Pi_+B\Pi_-\vert\psi\rangle,\qquad\langle\psi\vert\Pi_-B\Pi_+\vert\psi\rangle=\overline{z}
    (Π+BΠ−)†=Π−BΠ+(\Pi_+B\Pi_-)^\dagger=\Pi_-B\Pi_+ because BB and Π±\Pi_\pm are Hermitian, and ⟨ψ∣A†∣ψ⟩=⟨ψ∣A∣ψ⟩‾\langle\psi\vert A^\dagger\vert\psi\rangle=\overline{\langle\psi\vert A\vert\psi\rangle} for every operator AA.
  5. Step 5:
    f(θ)=a0+e2irθz+e−2irθz‾=a0+2 Re(e2irθz)f(\theta)=a_0+e^{2ir\theta}z+e^{-2ir\theta}\overline z=a_0+2\,\mathrm{Re}\big(e^{2ir\theta}z\big)
    Take the expectation value in ∣ψ⟩\lvert\psi\rangle; the two off-diagonal terms are complex conjugates, and w+w‾=2 Re ww+\overline w=2\,\mathrm{Re}\,w.
  6. Step 6:
    2 Re[(cos⁡2rθ+isin⁡2rθ)(Re z+i Im z)]=2 Re zcos⁡2rθ−2 Im zsin⁡2rθ2\,\mathrm{Re}\big[(\cos2r\theta+i\sin2r\theta)(\mathrm{Re}\,z+i\,\mathrm{Im}\,z)\big]=2\,\mathrm{Re}\,z\cos2r\theta-2\,\mathrm{Im}\,z\sin2r\theta
    Euler's formula and the real part of a product of two complex numbers, Re[(x+iy)(u+iv)]=xu−yv\mathrm{Re}[(x+iy)(u+iv)]=xu-yv.
  7. Step 7:
    a0=⟨ψ∣Π+BΠ+∣ψ⟩+⟨ψ∣Π−BΠ−∣ψ⟩∈Ra_0=\langle\psi\vert\Pi_+B\Pi_+\vert\psi\rangle+\langle\psi\vert\Pi_-B\Pi_-\vert\psi\rangle\in\mathbb R
    Π±BΠ±\Pi_\pm B\Pi_\pm is Hermitian, so its expectation value is real; a1a_1 and b1b_1 are real by construction.
  8. Step 8:
    G=12P:P 12(I±P)=12(P±I)=±12(I±P) ⇒ Π±=12(I±P), r=12G=\tfrac12P:\quad P\,\tfrac12(I\pm P)=\tfrac12(P\pm I)=\pm\tfrac12(I\pm P)\ \Rightarrow\ \Pi_\pm=\tfrac12(I\pm P),\ r=\tfrac12
    For a Pauli string P2=IP^2=I, so 12(I±P)\tfrac12(I\pm P) are the projectors onto the ±1\pm1 eigenspaces of PP, on which G=P/2G=P/2 takes the values ±12\pm\tfrac12. Then 2rθ=θ2r\theta=\theta.

∎

The shift identity and the general rule#

Theorem

General parameter-shift rule

Let f(θ)=a0+a1cos⁡2rθ+b1sin⁡2rθf(\theta)=a_0+a_1\cos2r\theta+b_1\sin2r\theta with r>0r>0. For every θ\theta and every shift ss, f(θ+s)−f(θ−s)=sin⁡2rsr f′(θ).f(\theta+s)-f(\theta-s)=\frac{\sin2rs}{r}\,f'(\theta). Consequently, whenever sin⁡2rs≠0\sin2rs\ne0, i.e. s∉π2rZs\notin\tfrac{\pi}{2r}\mathbb Z, f′(θ)=rsin⁡2rs[f(θ+s)−f(θ−s)].f'(\theta)=\frac{r}{\sin2rs}\Big[f(\theta+s)-f(\theta-s)\Big]. The shift s⋆=π4rs^\star=\frac{\pi}{4r} gives ∣sin⁡2rs∣=1\lvert\sin2rs\rvert=1, the smallest prefactor and hence the smallest amplification of statistical errors in the two evaluations; there f′(θ)=r[f(θ+π4r)−f(θ−π4r)]f'(\theta)=r\big[f(\theta+\tfrac{\pi}{4r})-f(\theta-\tfrac{\pi}{4r})\big].

Proof (7 steps)
  1. Step 1:
    f′(θ)=2r[−a1sin⁡2rθ+b1cos⁡2rθ]f'(\theta)=2r\big[-a_1\sin2r\theta+b_1\cos2r\theta\big]
    Differentiate term by term; the constant a0a_0 disappears and the chain rule produces the factor 2r2r.
  2. Step 2:
    cos⁡(x+y)−cos⁡(x−y)=−2sin⁡xsin⁡y,sin⁡(x+y)−sin⁡(x−y)=2cos⁡xsin⁡y\cos(x+y)-\cos(x-y)=-2\sin x\sin y,\qquad\sin(x+y)-\sin(x-y)=2\cos x\sin y
    Addition theorems: in cos⁡(x±y)=cos⁡xcos⁡y∓sin⁡xsin⁡y\cos(x\pm y)=\cos x\cos y\mp\sin x\sin y and sin⁡(x±y)=sin⁡xcos⁡y±cos⁡xsin⁡y\sin(x\pm y)=\sin x\cos y\pm\cos x\sin y the even parts cancel in the difference.
  3. Step 3:
    f(θ+s)−f(θ−s)=a1(−2sin⁡2rθ sin⁡2rs)+b1(2cos⁡2rθ sin⁡2rs)f(\theta+s)-f(\theta-s)=a_1\big(-2\sin2r\theta\,\sin2rs\big)+b_1\big(2\cos2r\theta\,\sin2rs\big)
    Apply both identities with x=2rθx=2r\theta and y=2rsy=2rs; the constant a0a_0 cancels in the difference.
  4. Step 4:
    =2sin⁡2rs [−a1sin⁡2rθ+b1cos⁡2rθ]=sin⁡2rsr f′(θ)=2\sin2rs\,\big[-a_1\sin2r\theta+b_1\cos2r\theta\big]=\frac{\sin2rs}{r}\,f'(\theta)
    Factor out 2sin⁡2rs2\sin2rs and compare with the first line: the bracket is f′(θ)/(2r)f'(\theta)/(2r).
  5. Step 5:
    sin⁡2rs≠0 ⇒ f′(θ)=rsin⁡2rs[f(θ+s)−f(θ−s)]\sin2rs\ne0\ \Rightarrow\ f'(\theta)=\frac{r}{\sin2rs}\big[f(\theta+s)-f(\theta-s)\big]
    Divide by sin⁡2rs/r\sin2rs/r. The identity holds exactly, for finite ss, with no remainder term.
  6. Step 6:
    Var[rsin⁡2rs(f^+−f^−)]=r2sin⁡22rs(σ+2+σ−2) ≥ r2(σ+2+σ−2)\mathrm{Var}\Big[\tfrac{r}{\sin2rs}\big(\hat f_+-\hat f_-\big)\Big]=\frac{r^2}{\sin^22rs}\big(\sigma_+^2+\sigma_-^2\big)\ \ge\ r^2\big(\sigma_+^2+\sigma_-^2\big)
    If the two values are estimated independently with variances σ±2\sigma_\pm^2, the prefactor enters squared. Equality holds exactly when ∣sin⁡2rs∣=1\lvert\sin2rs\rvert=1, i.e. s=π4r+kπ2rs=\tfrac{\pi}{4r}+\tfrac{k\pi}{2r}; the smallest positive choice is s⋆=π4rs^\star=\tfrac{\pi}{4r}.
  7. Step 7:
    s→0:rsin⁡2rs∼12s ⇒ f(θ+s)−f(θ−s)2s→f′(θ)s\to0:\quad\frac{r}{\sin2rs}\sim\frac{1}{2s}\ \Rightarrow\ \frac{f(\theta+s)-f(\theta-s)}{2s}\to f'(\theta)
    The central difference quotient is the small-shift limit of the rule; for finite ss it carries the factor sin⁡2rs2rs≠1\frac{\sin2rs}{2rs}\ne1, which is its truncation error. The shift rule replaces 12s\frac1{2s} by the exact factor rsin⁡2rs\frac{r}{\sin2rs}.

∎

The visualiser draws f(θ)=a0+a1cos⁡2rθ+b1sin⁡2rθf(\theta)=a_0+a_1\cos2r\theta+b_1\sin2r\theta, the two evaluation points θ±s\theta\pm s and the secant through them, next to the tangent at θ\theta. The secant slope is f(θ+s)−f(θ−s)2s=sin⁡2rs2rsf′(θ)\frac{f(\theta+s)-f(\theta-s)}{2s}=\frac{\sin2rs}{2rs}f'(\theta). It differs from the tangent slope by a fixed factor that does not depend on θ\theta or on the coefficients, and multiplying by 2rssin⁡2rs\frac{2rs}{\sin2rs} recovers the tangent exactly. The first instance shows a Pauli rotation (r=12r=\tfrac12, s=π2s=\tfrac\pi2); the second a generator with eigenvalues ±1\pm1, such as Z⊗ZZ\otimes Z in e−iθZ⊗Ze^{-i\theta Z\otimes Z}, with the small shift s=0.3s=0.3, where the secant is already close to the tangent and the correction factor is 0.6sin⁡0.6≈1.0626\frac{0.6}{\sin0.6}\approx1.0626.

Quick check

The gate e−iθ Z⊗Ze^{-i\theta\,Z\otimes Z} has the generator Z⊗ZZ\otimes Z with eigenvalues ±1\pm1. Which shift s∈(0,π2)s\in(0,\tfrac\pi2) minimises the variance of the shift-rule estimator rsin⁡2rs[f^(θ+s)−f^(θ−s)]\frac{r}{\sin2rs}\big[\hat f(\theta+s)-\hat f(\theta-s)\big]? Give ss in radians.

tolerance ±0.001

The two-term rule for Pauli rotations#

For a Pauli rotation RP(θ)=e−iθP/2R_P(\theta)=e^{-i\theta P/2} the generator P/2P/2 has eigenvalues ±12\pm\tfrac12, so r=12r=\tfrac12, the optimal shift is s⋆=π2s^\star=\tfrac\pi2 and the prefactor is rsin⁡2rs⋆=12\frac{r}{\sin2rs^\star}=\tfrac12. Applied to each parameter of a circuit in which every parameter occupies one Pauli rotation, this is the rule of Mitarai et al. (2018) and Schuld et al. (2019): ∂⟨H⟩∂θj=12[⟨H⟩θj+π2−⟨H⟩θj−π2].\frac{\partial\langle H\rangle}{\partial\theta_j}=\frac12\Big[\langle H\rangle_{\theta_j+\frac\pi2}-\langle H\rangle_{\theta_j-\frac\pi2}\Big]. Because this rule is the workhorse of variational quantum algorithms, we derive it twice more without the spectral detour. First we differentiate the state (Schrödinger picture), then the observable (Heisenberg picture). The second derivation yields an operator identity, which holds for every input state, mixed states included.

Derivation

Two-term rule in the Schrödinger picture

  1. Step 1:
    ∣ϕ(θ)⟩:=RP(θ)∣ψ⟩,f(θ)=⟨ϕ(θ)∣B∣ϕ(θ)⟩\lvert\phi(\theta)\rangle:=R_P(\theta)\lvert\psi\rangle,\qquad f(\theta)=\langle\phi(\theta)\vert B\vert\phi(\theta)\rangle
    In the Schrödinger picture the state carries the θ\theta-dependence and the observable BB is fixed.
  2. Step 2:
    ddθ∣ϕ(θ)⟩=−i2P e−iθP/2∣ψ⟩=−i2P ∣ϕ(θ)⟩\frac{d}{d\theta}\lvert\phi(\theta)\rangle=-\frac i2P\,e^{-i\theta P/2}\lvert\psi\rangle=-\frac i2P\,\lvert\phi(\theta)\rangle
    ddθe−iθG=−iGe−iθG\frac{d}{d\theta}e^{-i\theta G}=-iGe^{-i\theta G} for a fixed operator GG, here G=P/2G=P/2.
  3. Step 3:
    f′(θ)=(i2⟨ϕ∣P)B∣ϕ⟩+⟨ϕ∣B(−i2P∣ϕ⟩)=i2⟨ϕ∣[P,B]∣ϕ⟩f'(\theta)=\Big(\frac i2\langle\phi\vert P\Big)B\vert\phi\rangle+\langle\phi\vert B\Big(-\frac i2P\vert\phi\rangle\Big)=\frac i2\langle\phi\vert[P,B]\vert\phi\rangle
    Product rule for the bra and the ket. The bra of −i2P∣ϕ⟩-\tfrac i2P\lvert\phi\rangle is i2⟨ϕ∣P\tfrac i2\langle\phi\rvert P, because PP is Hermitian and complex conjugation turns −i-i into ii.
  4. Step 4:
    RP(±π2)=cos⁡π4 I∓isin⁡π4 P=12(I∓iP)R_P\big(\pm\tfrac\pi2\big)=\cos\tfrac\pi4\,I\mp i\sin\tfrac\pi4\,P=\tfrac1{\sqrt2}\big(I\mp iP\big)
    P2=IP^2=I gives e−iαP/2=cos⁡α2 I−isin⁡α2 Pe^{-i\alpha P/2}=\cos\frac\alpha2\,I-i\sin\frac\alpha2\,P (Euler's formula for an involution), evaluated at α=±π2\alpha=\pm\frac\pi2.
  5. Step 5:
    RP(θ±π2)∣ψ⟩=RP(±π2)RP(θ)∣ψ⟩=12(I∓iP)∣ϕ(θ)⟩R_P\big(\theta\pm\tfrac\pi2\big)\lvert\psi\rangle=R_P\big(\pm\tfrac\pi2\big)R_P(\theta)\lvert\psi\rangle=\tfrac1{\sqrt2}\big(I\mp iP\big)\lvert\phi(\theta)\rangle
    Rotations about the same axis compose additively: e−iaP/2e−ibP/2=e−i(a+b)P/2e^{-iaP/2}e^{-ibP/2}=e^{-i(a+b)P/2}, because the exponents commute.
  6. Step 6:
    f(θ±π2)=12⟨ϕ∣(I±iP) B (I∓iP)∣ϕ⟩=12[⟨B⟩+⟨PBP⟩±i⟨PB⟩∓i⟨BP⟩]f\big(\theta\pm\tfrac\pi2\big)=\tfrac12\langle\phi\vert(I\pm iP)\,B\,(I\mp iP)\vert\phi\rangle=\tfrac12\Big[\langle B\rangle+\langle PBP\rangle\pm i\langle PB\rangle\mp i\langle BP\rangle\Big]
    The bra of (I∓iP)∣ϕ⟩(I\mp iP)\lvert\phi\rangle is ⟨ϕ∣(I±iP)\langle\phi\rvert(I\pm iP); expand the product, with ⟨⋅⟩\langle\cdot\rangle the expectation value in ∣ϕ(θ)⟩\lvert\phi(\theta)\rangle.
  7. Step 7:
    f(θ+π2)−f(θ−π2)=i⟨PB⟩−i⟨BP⟩=i⟨ϕ∣[P,B]∣ϕ⟩f\big(\theta+\tfrac\pi2\big)-f\big(\theta-\tfrac\pi2\big)=i\langle PB\rangle-i\langle BP\rangle=i\langle\phi\vert[P,B]\vert\phi\rangle
    The terms ⟨B⟩\langle B\rangle and ⟨PBP⟩\langle PBP\rangle are common to both shifts and cancel; the cross terms appear with opposite signs and add.
  8. Step 8:
    f′(θ)=i2⟨[P,B]⟩=12[f(θ+π2)−f(θ−π2)]f'(\theta)=\frac i2\langle[P,B]\rangle=\frac12\Big[f\big(\theta+\tfrac\pi2\big)-f\big(\theta-\tfrac\pi2\big)\Big]
    Compare the third and the seventh line: both are multiples of the same commutator expectation value.

Result

ddθ⟨ψ∣RP†(θ)BRP(θ)∣ψ⟩=12[f(θ+π2)−f(θ−π2)]\boxed{\frac{d}{d\theta}\langle\psi\vert R_P^\dagger(\theta)BR_P(\theta)\vert\psi\rangle=\frac12\Big[f\big(\theta+\tfrac\pi2\big)-f\big(\theta-\tfrac\pi2\big)\Big]}
Proposition

Two-term rule in the Heisenberg picture

Let U(θ)=e−iθGU(\theta)=e^{-i\theta G} with G=12PG=\tfrac12P, PP a Pauli string, and let BB be Hermitian. Then ddθ[U†(θ)BU(θ)]=iU†(θ)[G,B]U(θ)=12[U†(θ+π2)BU(θ+π2)−U†(θ−π2)BU(θ−π2)].\frac{d}{d\theta}\Big[U^\dagger(\theta)BU(\theta)\Big]=iU^\dagger(\theta)[G,B]U(\theta)=\frac12\Big[U^\dagger\big(\theta+\tfrac\pi2\big)BU\big(\theta+\tfrac\pi2\big)-U^\dagger\big(\theta-\tfrac\pi2\big)BU\big(\theta-\tfrac\pi2\big)\Big]. The commutator is itself a difference of two conjugations, i[G,B]=12[U†(π2)BU(π2)−U†(−π2)BU(−π2)]i[G,B]=\tfrac12\big[U^\dagger(\tfrac\pi2)BU(\tfrac\pi2)-U^\dagger(-\tfrac\pi2)BU(-\tfrac\pi2)\big]. Taking the expectation value in any state, pure or mixed, gives the two-term rule.

Proof (8 steps)
  1. Step 1:
    B(θ):=U†(θ)BU(θ),f(θ)=⟨ψ∣B(θ)∣ψ⟩B(\theta):=U^\dagger(\theta)BU(\theta),\qquad f(\theta)=\langle\psi\vert B(\theta)\vert\psi\rangle
    In the Heisenberg picture the observable carries the θ\theta-dependence and the state is fixed.
  2. Step 2:
    dUdθ=−iGU=−iUG,dU†dθ=iU†G=iGU†\frac{dU}{d\theta}=-iGU=-iUG,\qquad\frac{dU^\dagger}{d\theta}=iU^\dagger G=iGU^\dagger
    GG commutes with every function of itself, in particular with e∓iθGe^{\mp i\theta G}.
  3. Step 3:
    dB(θ)dθ=iU†GBU−iU†BGU=iU†[G,B]U\frac{dB(\theta)}{d\theta}=iU^\dagger GBU-iU^\dagger BGU=iU^\dagger[G,B]U
    Product rule for the three factors U†U^\dagger, BB, UU; the middle factor is constant.
  4. Step 4:
    U(±π2)=12(I∓iP),U†(±π2)=12(I±iP)U\big(\pm\tfrac\pi2\big)=\tfrac1{\sqrt2}(I\mp iP),\qquad U^\dagger\big(\pm\tfrac\pi2\big)=\tfrac1{\sqrt2}(I\pm iP)
    The same Euler formula as in the Schrödinger derivation, with P†=PP^\dagger=P for the adjoint.
  5. Step 5:
    U†(±π2) B U(±π2)=12(B+PBP±iPB∓iBP)U^\dagger\big(\pm\tfrac\pi2\big)\,B\,U\big(\pm\tfrac\pi2\big)=\tfrac12\big(B+PBP\pm iPB\mp iBP\big)
    Expand 12(I±iP)B(I∓iP)=12(B∓iBP±iPB+PBP)\tfrac12(I\pm iP)B(I\mp iP)=\tfrac12\big(B\mp iBP\pm iPB+PBP\big), using (±i)(∓i)=1(\pm i)(\mp i)=1.
  6. Step 6:
    U†(π2)BU(π2)−U†(−π2)BU(−π2)=i(PB−BP)=i[P,B]=2i[G,B]U^\dagger\big(\tfrac\pi2\big)BU\big(\tfrac\pi2\big)-U^\dagger\big(-\tfrac\pi2\big)BU\big(-\tfrac\pi2\big)=i(PB-BP)=i[P,B]=2i[G,B]
    Subtract the two signs of the previous line; B+PBPB+PBP cancels, and [P,B]=2[G,B][P,B]=2[G,B] because G=P/2G=P/2.
  7. Step 7:
    U†(θ) U†(±π2) B U(±π2) U(θ)=U†(θ±π2) B U(θ±π2)U^\dagger(\theta)\,U^\dagger\big(\pm\tfrac\pi2\big)\,B\,U\big(\pm\tfrac\pi2\big)\,U(\theta)=U^\dagger\big(\theta\pm\tfrac\pi2\big)\,B\,U\big(\theta\pm\tfrac\pi2\big)
    U(a)U(b)=U(a+b)U(a)U(b)=U(a+b) for the one-parameter group e−iθGe^{-i\theta G}, and the adjoint reverses the order.
  8. Step 8:
    f′(θ)=⟨ψ∣iU†[G,B]U∣ψ⟩=12[f(θ+π2)−f(θ−π2)],ddθTr[ρB(θ)]=12[Tr[ρB(θ+π2)]−Tr[ρB(θ−π2)]]f'(\theta)=\langle\psi\vert iU^\dagger[G,B]U\vert\psi\rangle=\frac12\Big[f\big(\theta+\tfrac\pi2\big)-f\big(\theta-\tfrac\pi2\big)\Big],\qquad\frac{d}{d\theta}\mathrm{Tr}\big[\rho B(\theta)\big]=\frac12\Big[\mathrm{Tr}\big[\rho B(\theta+\tfrac\pi2)\big]-\mathrm{Tr}\big[\rho B(\theta-\tfrac\pi2)\big]\Big]
    Substitute line six into line three, conjugate by U(θ)U(\theta) with line seven and take the expectation value. The identity holds between operators, so it holds for every density matrix ρ\rho as well.

∎

The circuit of the lab's visible test makes the rule concrete. Two RYR_Y rotations prepare a product state and a CNOT correlates the qubits before Z1Z_1 is measured.

Circuit · 2 qubits

Quantum circuit on 2 qubitsRY(0.40) on q0; RY(1.10) on q1; X on q1 controlled by q0q0q1RY(0.40)RY(1.10)

The lab's visible circuit at (t1,t2)=(0.4,1.1)(t_1,t_2)=(0.4,1.1): RY(t1)R_Y(t_1) on qubit 00, RY(t2)R_Y(t_2) on qubit 11, then CNOT(0,1)\mathrm{CNOT}(0,1). The observable is Z1Z_1 (Pauli word IZ\texttt{IZ}), and ⟨Z1⟩=cos⁡0.4cos⁡1.1≈0.4178\langle Z_1\rangle=\cos0.4\cos1.1\approx0.4178.

∣ψ⟩=0.8355∣00⟩+0.5123∣01⟩+0.1038∣10⟩+0.1694∣11⟩\lvert\psi\rangle = 0.8355\lvert 00\rangle + 0.5123\lvert 01\rangle + 0.1038\lvert 10\rangle + 0.1694\lvert 11\rangle
Derivation

The visible lab test: f = cos t1 cos t2

  1. Step 1:
    RY(t1)⊗RY(t2)∣00⟩=(c1∣0⟩+s1∣1⟩)⊗(c2∣0⟩+s2∣1⟩),ck=cos⁡tk2, sk=sin⁡tk2R_Y(t_1)\otimes R_Y(t_2)\lvert00\rangle=\big(c_1\lvert0\rangle+s_1\lvert1\rangle\big)\otimes\big(c_2\lvert0\rangle+s_2\lvert1\rangle\big),\qquad c_k=\cos\tfrac{t_k}2,\ s_k=\sin\tfrac{t_k}2
    The first column of RY(t)=(cos⁡t2−sin⁡t2sin⁡t2cos⁡t2)R_Y(t)=\begin{pmatrix}\cos\frac t2&-\sin\frac t2\\\sin\frac t2&\cos\frac t2\end{pmatrix} on each qubit.
  2. Step 2:
    CNOT† Z1 CNOT=Z0Z1\mathrm{CNOT}^\dagger\,Z_1\,\mathrm{CNOT}=Z_0Z_1
    CNOT∣a,b⟩=∣a,a⊕b⟩\mathrm{CNOT}\lvert a,b\rangle=\lvert a,a\oplus b\rangle, and Z1Z_1 then has the eigenvalue (−1)a⊕b=(−1)a(−1)b(-1)^{a\oplus b}=(-1)^a(-1)^b, which is the eigenvalue of Z0Z1Z_0Z_1 on ∣a,b⟩\lvert a,b\rangle.
  3. Step 3:
    f(t1,t2)=⟨Z⟩q0 ⟨Z⟩q1=(c12−s12)(c22−s22)=cos⁡t1cos⁡t2f(t_1,t_2)=\langle Z\rangle_{q_0}\,\langle Z\rangle_{q_1}=\big(c_1^2-s_1^2\big)\big(c_2^2-s_2^2\big)=\cos t_1\cos t_2
    On a product state the expectation value of a product of local operators factorises; cos⁡2t2−sin⁡2t2=cos⁡t\cos^2\frac t2-\sin^2\frac t2=\cos t.
  4. Step 4:
    f(t1±π2,t2)=cos⁡(t1±π2)cos⁡t2=∓sin⁡t1cos⁡t2f\big(t_1\pm\tfrac\pi2,t_2\big)=\cos\big(t_1\pm\tfrac\pi2\big)\cos t_2=\mp\sin t_1\cos t_2
    cos⁡(x±π2)=∓sin⁡x\cos(x\pm\frac\pi2)=\mp\sin x.
  5. Step 5:
    12[f(t1+π2,t2)−f(t1−π2,t2)]=12(−sin⁡t1cos⁡t2−sin⁡t1cos⁡t2)=−sin⁡t1cos⁡t2=∂f∂t1\tfrac12\Big[f\big(t_1+\tfrac\pi2,t_2\big)-f\big(t_1-\tfrac\pi2,t_2\big)\Big]=\tfrac12\big(-\sin t_1\cos t_2-\sin t_1\cos t_2\big)=-\sin t_1\cos t_2=\frac{\partial f}{\partial t_1}
    The two-term rule reproduces the analytic partial derivative; the same computation in t2t_2 gives −cos⁡t1sin⁡t2-\cos t_1\sin t_2.
  6. Step 6:
    ∇f(0.4, 1.1)=(−sin⁡0.4cos⁡1.1, −cos⁡0.4sin⁡1.1)=(−0.1766386497, −0.8208563369)\nabla f(0.4,\,1.1)=\big(-\sin0.4\cos1.1,\ -\cos0.4\sin1.1\big)=(-0.1766386497,\ -0.8208563369)
    The expected value of the visible test, obtained with four evaluations, two per parameter.

Quick check

The parameter θ3\theta_3 of a circuit sits in a single RXR_X gate. Two runs give f(θ+π2e3)=0.42f(\theta+\tfrac\pi2e_3)=0.42 and f(θ−π2e3)=−0.18f(\theta-\tfrac\pi2e_3)=-0.18. What is ∂f/∂θ3\partial f/\partial\theta_3 at θ\theta?

tolerance ±0.000001

Exact, and unbiased under shot noise#

On hardware, ff is not available exactly. Each evaluation is an average over NN measurement shots. Two kinds of error must be separated: the systematic error of the formula, which remains even with infinitely many shots, and the statistical error of the shot averages. Finite differences suffer from both. The parameter-shift rule has no systematic error at all.

Theorem

Bias and variance of the shift and difference estimators

Let ff be the cost of a Pauli rotation (r=12r=\tfrac12), measured with a Pauli observable BB (eigenvalues ±1\pm1), and let f^(ϑ)\hat f(\vartheta) be the mean of NN independent single-shot outcomes at angle ϑ\vartheta. (a) The central difference quotient with step hh satisfies f(θ+h)−f(θ−h)2h=sin⁡hhf′(θ)\frac{f(\theta+h)-f(\theta-h)}{2h}=\frac{\sin h}{h}f'(\theta), so its bias is (sin⁡hh−1)f′(θ)=−h26f′(θ)+O(h4)\big(\frac{\sin h}{h}-1\big)f'(\theta)=-\frac{h^2}{6}f'(\theta)+O(h^4). (b) The parameter-shift estimator g^=12[f^(θ+π2)−f^(θ−π2)]\hat g=\tfrac12\big[\hat f(\theta+\tfrac\pi2)-\hat f(\theta-\tfrac\pi2)\big] is unbiased, E[g^]=f′(θ)\mathbb E[\hat g]=f'(\theta), with Var[g^]=(1−f+2)+(1−f−2)4N≤12N,f±=f(θ±π2).\mathrm{Var}[\hat g]=\frac{\big(1-f_+^2\big)+\big(1-f_-^2\big)}{4N}\le\frac1{2N},\qquad f_\pm=f\big(\theta\pm\tfrac\pi2\big). (c) The central-difference estimator g^h=f^(θ+h)−f^(θ−h)2h\hat g_h=\frac{\hat f(\theta+h)-\hat f(\theta-h)}{2h} has Var[g^h]=(1−f(θ+h)2)+(1−f(θ−h)2)4h2N\mathrm{Var}[\hat g_h]=\frac{(1-f(\theta+h)^2)+(1-f(\theta-h)^2)}{4h^2N}.

Proof (8 steps)
  1. Step 1:
    f(θ+h)−f(θ−h)=sin⁡(2⋅12h)12f′(θ)=2sin⁡h f′(θ)f(\theta+h)-f(\theta-h)=\frac{\sin(2\cdot\frac12h)}{\frac12}f'(\theta)=2\sin h\,f'(\theta)
    The shift identity of the general rule with r=12r=\tfrac12 and s=hs=h.
  2. Step 2:
    f(θ+h)−f(θ−h)2h=sin⁡hhf′(θ),sin⁡hh=1−h26+O(h4)\frac{f(\theta+h)-f(\theta-h)}{2h}=\frac{\sin h}{h}f'(\theta),\qquad\frac{\sin h}{h}=1-\frac{h^2}{6}+O(h^4)
    Divide by 2h2h and expand sin⁡h=h−h36+O(h5)\sin h=h-\frac{h^3}6+O(h^5). The truncation error of finite differences is therefore an exact multiplicative factor for such circuits, and it vanishes only as h→0h\to0.
  3. Step 3:
    X∈{+1,−1},Pr⁡(X=±1)=12(1±f),E[X]=f,Var[X]=1−f2X\in\{+1,-1\},\quad\Pr(X=\pm1)=\tfrac12(1\pm f),\quad\mathbb E[X]=f,\quad\mathrm{Var}[X]=1-f^2
    Born rule for an observable with eigenvalues ±1\pm1: the probabilities are ⟨12(I±B)⟩\langle\tfrac12(I\pm B)\rangle. Since X2=1X^2=1, Var[X]=1−f2\mathrm{Var}[X]=1-f^2.
  4. Step 4:
    f^=1N∑k=1NXk:E[f^]=f,Var[f^]=1−f2N\hat f=\frac1N\sum_{k=1}^NX_k:\qquad\mathbb E[\hat f]=f,\qquad\mathrm{Var}[\hat f]=\frac{1-f^2}{N}
    Linearity of expectation and additivity of variance for independent, identically distributed shots.
  5. Step 5:
    E[g^]=12(f+−f−)=f′(θ)\mathbb E[\hat g]=\tfrac12\big(f_+-f_-\big)=f'(\theta)
    Linearity of expectation, followed by the two-term rule, which holds exactly. No step size appears, so no truncation error survives in the expectation.
  6. Step 6:
    Var[g^]=14(Var[f^(θ+π2)]+Var[f^(θ−π2)])=(1−f+2)+(1−f−2)4N≤24N\mathrm{Var}[\hat g]=\tfrac14\Big(\mathrm{Var}\big[\hat f(\theta+\tfrac\pi2)\big]+\mathrm{Var}\big[\hat f(\theta-\tfrac\pi2)\big]\Big)=\frac{(1-f_+^2)+(1-f_-^2)}{4N}\le\frac{2}{4N}
    The two circuit runs are independent, so their variances add, each scaled by (12)2(\tfrac12)^2. The bound uses f±2≥0f_\pm^2\ge0.
  7. Step 7:
    Var[g^h]=14h2(Var[f^(θ+h)]+Var[f^(θ−h)])=(1−f(θ+h)2)+(1−f(θ−h)2)4h2N\mathrm{Var}[\hat g_h]=\frac{1}{4h^2}\Big(\mathrm{Var}\big[\hat f(\theta+h)\big]+\mathrm{Var}\big[\hat f(\theta-h)\big]\Big)=\frac{(1-f(\theta+h)^2)+(1-f(\theta-h)^2)}{4h^2N}
    The same computation with the prefactor 12h\frac1{2h} in place of 12\tfrac12: the variance grows like h−2h^{-2}, which is why a small step makes finite differences useless under shot noise.
  8. Step 8:
    f=cos⁡θ, θ=π4, N=1000:Var[g^]=12+124000=2.5×10−4,Var[g^0.01]≈14⋅10−4⋅1000=2.5f=\cos\theta,\ \theta=\tfrac\pi4,\ N=1000:\qquad\mathrm{Var}[\hat g]=\frac{\frac12+\frac12}{4000}=2.5\times10^{-4},\qquad\mathrm{Var}[\hat g_{0.01}]\approx\frac{1}{4\cdot10^{-4}\cdot1000}=2.5
    The example of the previous lesson: f±=cos⁡(π4±π2)=∓22f_\pm=\cos(\frac\pi4\pm\frac\pi2)=\mp\frac{\sqrt2}2, while f(π4±h)2≈12f(\frac\pi4\pm h)^2\approx\frac12. With the same shot budget the shift estimator is 10410^4 times less noisy than central differences with h=0.01h=0.01.

∎

Quick check

A Pauli observable is measured with N=1000N=1000 shots at each of the two shifted points, where its exact values are f(θ+π2)=0.6f(\theta+\tfrac\pi2)=0.6 and f(θ−π2)=−0.6f(\theta-\tfrac\pi2)=-0.6. What is the variance of the parameter-shift estimator g^\hat g?

tolerance ±0.000001

Parameters that occur in several gates#

In data re-uploading circuits, in QAOA and in ansätze with tied weights, one trainable parameter α\alpha feeds several gates, often through an affine map such as RZZ(2γ)R_{ZZ}(2\gamma). The shift rule applies to gate angles, not directly to parameters. The chain rule connects the two.

Proposition

Chain rule for repeated parameters

Let a parameter α\alpha enter Pauli rotations g=1,…,mg=1,\dots,m with angles θg=agα+bg\theta_g=a_g\alpha+b_g, and let f(θ1,…,θm)f(\theta_1,\dots,\theta_m) be the expectation value as a function of independent gate angles, so that F(α)=f(θ1(α),…,θm(α))F(\alpha)=f(\theta_1(\alpha),\dots,\theta_m(\alpha)). Then dFdα=∑g=1mag⋅12[f(θ+π2eg)−f(θ−π2eg)],\frac{dF}{d\alpha}=\sum_{g=1}^ma_g\cdot\frac12\Big[f\big(\theta+\tfrac\pi2e_g\big)-f\big(\theta-\tfrac\pi2e_g\big)\Big], where θ±π2eg\theta\pm\tfrac\pi2e_g shifts the angle of occurrence gg only. The gradient costs 2m2m evaluations, and occurrences with ag=0a_g=0 can be skipped.

Proof (6 steps)
  1. Step 1:
    F(α)=f(θ1(α),…,θm(α)),θg(α)=agα+bgF(\alpha)=f\big(\theta_1(\alpha),\dots,\theta_m(\alpha)\big),\qquad\theta_g(\alpha)=a_g\alpha+b_g
    Give each occurrence its own angle variable; FF is ff restricted to the line θ(α)\theta(\alpha) in angle space.
  2. Step 2:
    dFdα=∑g=1m∂f∂θg dθgdα=∑g=1mag ∂f∂θg\frac{dF}{d\alpha}=\sum_{g=1}^m\frac{\partial f}{\partial\theta_g}\,\frac{d\theta_g}{d\alpha}=\sum_{g=1}^ma_g\,\frac{\partial f}{\partial\theta_g}
    Multivariate chain rule; ff is smooth (a trigonometric polynomial) in every angle.
  3. Step 3:
    θg↦f(…,θg,… )=a0(g)+a1(g)cos⁡θg+b1(g)sin⁡θg\theta_g\mapsto f(\dots,\theta_g,\dots)=a_0^{(g)}+a_1^{(g)}\cos\theta_g+b_1^{(g)}\sin\theta_g
    With all other angles fixed, the remaining gates are absorbed into ∣ψ⟩\lvert\psi\rangle and BB as in the first derivation, and the single-frequency theorem applies to occurrence gg alone.
  4. Step 4:
    ∂f∂θg=12[f(θ+π2eg)−f(θ−π2eg)]\frac{\partial f}{\partial\theta_g}=\frac12\Big[f\big(\theta+\tfrac\pi2e_g\big)-f\big(\theta-\tfrac\pi2e_g\big)\Big]
    The two-term rule for occurrence gg. The other occurrences of α\alpha keep their unshifted angles.
  5. Step 5:
    dFdα=∑gag 12[f(θ+π2eg)−f(θ−π2eg)],#evaluations=2 #{g:ag≠0}\frac{dF}{d\alpha}=\sum_ga_g\,\frac12\Big[f\big(\theta+\tfrac\pi2e_g\big)-f\big(\theta-\tfrac\pi2e_g\big)\Big],\qquad\#\text{evaluations}=2\,\#\{g:a_g\ne0\}
    Insert line four into line two. Shifting the gate angle by ±π2\pm\tfrac\pi2 is the same as shifting α\alpha by ±π2ag\pm\frac{\pi}{2a_g} in that one gate.
  6. Step 6:
    RZZ(2γ): a=2 ⇒ ∂F∂γ=2⋅12[f(θZZ+π2)−f(θZZ−π2)]R_{ZZ}(2\gamma):\ a=2\ \Rightarrow\ \frac{\partial F}{\partial\gamma}=2\cdot\frac12\Big[f\big(\theta_{ZZ}+\tfrac\pi2\big)-f\big(\theta_{ZZ}-\tfrac\pi2\big)\Big]
    A single occurrence with scale a=2a=2: the derivative with respect to the parameter is twice the derivative with respect to the gate angle.

∎

Derivation

Shifting all occurrences at once is wrong

  1. Step 1:
    RY(α)RY(α)∣0⟩=RY(2α)∣0⟩,F(α)=⟨Z⟩=cos⁡2α,F′(α)=−2sin⁡2αR_Y(\alpha)R_Y(\alpha)\lvert0\rangle=R_Y(2\alpha)\lvert0\rangle,\qquad F(\alpha)=\langle Z\rangle=\cos2\alpha,\qquad F'(\alpha)=-2\sin2\alpha
    Two rotations about the same axis add their angles, and ⟨Z⟩\langle Z\rangle of RY(β)∣0⟩R_Y(\beta)\lvert0\rangle is cos⁡β\cos\beta.
  2. Step 2:
    12[F(α+π2)−F(α−π2)]=12[cos⁡(2α+π)−cos⁡(2α−π)]=12[−cos⁡2α+cos⁡2α]=0\tfrac12\Big[F\big(\alpha+\tfrac\pi2\big)-F\big(\alpha-\tfrac\pi2\big)\Big]=\tfrac12\big[\cos(2\alpha+\pi)-\cos(2\alpha-\pi)\big]=\tfrac12\big[-\cos2\alpha+\cos2\alpha\big]=0
    Shifting α\alpha shifts both gate angles, i.e. the total angle by ±π\pm\pi, and cos⁡(x±π)=−cos⁡x\cos(x\pm\pi)=-\cos x. The result is 00 for every α\alpha, which is wrong.
  3. Step 3:
    f(θ1,θ2)=cos⁡(θ1+θ2),∂f∂θg=12[cos⁡(θ1+θ2+π2)−cos⁡(θ1+θ2−π2)]=−sin⁡(θ1+θ2)f(\theta_1,\theta_2)=\cos(\theta_1+\theta_2),\qquad\frac{\partial f}{\partial\theta_g}=\tfrac12\Big[\cos\big(\theta_1+\theta_2+\tfrac\pi2\big)-\cos\big(\theta_1+\theta_2-\tfrac\pi2\big)\Big]=-\sin(\theta_1+\theta_2)
    With independent occurrence angles, each shift moves the total angle by only ±π2\pm\tfrac\pi2; cos⁡(x±π2)=∓sin⁡x\cos(x\pm\frac\pi2)=\mp\sin x.
  4. Step 4:
    dFdα=1⋅(−sin⁡2α)+1⋅(−sin⁡2α)=−2sin⁡2α\frac{dF}{d\alpha}=1\cdot\big(-\sin2\alpha\big)+1\cdot\big(-\sin2\alpha\big)=-2\sin2\alpha
    The chain rule with a1=a2=1a_1=a_2=1, evaluated at θ1=θ2=α\theta_1=\theta_2=\alpha, reproduces F′(α)F'(\alpha) with four evaluations.
  5. Step 5:
    F(α)=cos⁡2α: r=1, s⋆=π4 ⇒ F(α+π4)−F(α−π4)=−sin⁡2α−sin⁡2α=−2sin⁡2αF(\alpha)=\cos2\alpha:\ r=1,\ s^\star=\tfrac\pi4\ \Rightarrow\ F\big(\alpha+\tfrac\pi4\big)-F\big(\alpha-\tfrac\pi4\big)=-\sin2\alpha-\sin2\alpha=-2\sin2\alpha
    Alternatively, as a function of α\alpha the circuit is e−iαYe^{-i\alpha Y}, a generator with eigenvalues ±1\pm1, so the general rule with r=1r=1 and shift π4\tfrac\pi4 applies (two evaluations). The naive rule used r=12r=\tfrac12 and s=π2s=\tfrac\pi2, where sin⁡2rs=sin⁡π=0\sin2rs=\sin\pi=0.

Quick check

A parameter α\alpha appears in RX(α)R_X(\alpha) on qubit 00 and in RZ(3α)R_Z(3\alpha) on qubit 11 of the same circuit. Which procedure returns dF/dαdF/d\alpha exactly?

Select one answer

Controlled rotations and the four-term rule#

A controlled rotation CRP(θ)CR_P(\theta) applies RP(θ)R_P(\theta) to the target when the control is ∣1⟩\lvert1\rangle and does nothing otherwise. It is a one-parameter gate, but its generator has three distinct eigenvalues. The two-term rule then fails, and the correct rule needs four evaluations (Anselmetti et al. 2021; Wierichs et al. 2022).

Theorem

Spectrum and frequencies of a controlled rotation

Let CRP(θ)=∣0⟩⟨0∣⊗I+∣1⟩⟨1∣⊗RP(θ)CR_P(\theta)=\lvert0\rangle\langle0\rvert\otimes I+\lvert1\rangle\langle1\rvert\otimes R_P(\theta) with control qubit first. Then CRP(θ)=e−iθGCR_P(\theta)=e^{-i\theta G} with G=∣1⟩⟨1∣⊗12PG=\lvert1\rangle\langle1\rvert\otimes\tfrac12P, whose eigenvalues are 00 (twice), 12\tfrac12 and −12-\tfrac12, and f(θ)=a0+a1cos⁡θ2+b1sin⁡θ2+a2cos⁡θ+b2sin⁡θf(\theta)=a_0+a_1\cos\tfrac\theta2+b_1\sin\tfrac\theta2+a_2\cos\theta+b_2\sin\theta with real coefficients. With Π0=∣0⟩⟨0∣⊗I\Pi_0=\lvert0\rangle\langle0\rvert\otimes I, Π±=∣1⟩⟨1∣⊗12(I±P)\Pi_\pm=\lvert1\rangle\langle1\rvert\otimes\tfrac12(I\pm P) and zjk=⟨ψ∣ΠjBΠk∣ψ⟩z_{jk}=\langle\psi\vert\Pi_jB\Pi_k\vert\psi\rangle: a0=z00+z+++z−−a_0=z_{00}+z_{++}+z_{--}, a1−ib1=2(z+0+z0−)a_1-ib_1=2(z_{+0}+z_{0-}) and a2−ib2=2z+−a_2-ib_2=2z_{+-}.

Proof (7 steps)
  1. Step 1:
    G=∣1⟩⟨1∣⊗12P=(00012P) ⇒ e−iθG=(I00e−iθP/2)=CRP(θ)G=\lvert1\rangle\langle1\rvert\otimes\tfrac12P=\begin{pmatrix}0&0\\0&\tfrac12P\end{pmatrix}\ \Rightarrow\ e^{-i\theta G}=\begin{pmatrix}I&0\\0&e^{-i\theta P/2}\end{pmatrix}=CR_P(\theta)
    GG is block diagonal in the control basis, and the exponential of a block-diagonal matrix is taken block by block: e0=Ie^{0}=I on the control-00 block.
  2. Step 2:
    G ∣0⟩∣v⟩=0,G ∣1⟩∣p±⟩=±12 ∣1⟩∣p±⟩,P∣p±⟩=±∣p±⟩G\,\lvert0\rangle\lvert v\rangle=0,\qquad G\,\lvert1\rangle\lvert p_\pm\rangle=\pm\tfrac12\,\lvert1\rangle\lvert p_\pm\rangle,\qquad P\lvert p_\pm\rangle=\pm\lvert p_\pm\rangle
    The control-00 block contributes the eigenvalue 00 with multiplicity 22; the control-11 block contributes the eigenvalues ±12\pm\tfrac12 of 12P\tfrac12P.
  3. Step 3:
    G=0⋅Π0+12Π+−12Π−,U(θ)=Π0+e−iθ/2Π++eiθ/2Π−G=0\cdot\Pi_0+\tfrac12\Pi_+-\tfrac12\Pi_-,\qquad U(\theta)=\Pi_0+e^{-i\theta/2}\Pi_++e^{i\theta/2}\Pi_-
    Spectral decomposition with three projectors summing to II, and the operator function e−iθλe^{-i\theta\lambda} applied eigenvalue by eigenvalue.
  4. Step 4:
    U†BU=∑j,k∈{0,+,−}ei(λj−λk)θ ΠjBΠk,λ0=0, λ±=±12U^\dagger BU=\sum_{j,k\in\{0,+,-\}}e^{i(\lambda_j-\lambda_k)\theta}\,\Pi_jB\Pi_k,\qquad\lambda_0=0,\ \lambda_\pm=\pm\tfrac12
    Multiply out U†=∑jeiλjθΠjU^\dagger=\sum_je^{i\lambda_j\theta}\Pi_j and U=∑ke−iλkθΠkU=\sum_ke^{-i\lambda_k\theta}\Pi_k, as in the two-eigenvalue case.
  5. Step 5:
    λj−λk∈{0,±12,±1}:12 from (+,0),(0,−);1 from (+,−);0 from j=k\lambda_j-\lambda_k\in\{0,\pm\tfrac12,\pm1\}:\quad\tfrac12\text{ from }(+,0),(0,-);\qquad1\text{ from }(+,-);\qquad0\text{ from }j=k
    The frequencies of ff are the differences of generator eigenvalues. Two different pairs produce the frequency 12\tfrac12, and −12-\tfrac12 arises from the reversed pairs.
  6. Step 6:
    f(θ)=a0+2 Re[eiθ/2(z+0+z0−)]+2 Re[eiθz+−],zkj=zjk‾f(\theta)=a_0+2\,\mathrm{Re}\Big[e^{i\theta/2}\big(z_{+0}+z_{0-}\big)\Big]+2\,\mathrm{Re}\Big[e^{i\theta}z_{+-}\Big],\qquad z_{kj}=\overline{z_{jk}}
    Group each term with its reversed partner, which is its complex conjugate because (ΠjBΠk)†=ΠkBΠj(\Pi_jB\Pi_k)^\dagger=\Pi_kB\Pi_j.
  7. Step 7:
    2 Re(eiωθw)=2 Re wcos⁡ωθ−2 Im wsin⁡ωθ2\,\mathrm{Re}\big(e^{i\omega\theta}w\big)=2\,\mathrm{Re}\,w\cos\omega\theta-2\,\mathrm{Im}\,w\sin\omega\theta
    As in the two-eigenvalue theorem, with ω=12\omega=\tfrac12, w=z+0+z0−w=z_{+0}+z_{0-} and with ω=1\omega=1, w=z+−w=z_{+-}; this gives the stated a1,b1,a2,b2a_1,b_1,a_2,b_2.

∎

Derivation

Why the two-term rule is wrong for controlled rotations

  1. Step 1:
    f=a0+f1/2+f1,fω(θ)=aωcos⁡ωθ+bωsin⁡ωθf=a_0+f_{1/2}+f_1,\qquad f_\omega(\theta)=a_\omega\cos\omega\theta+b_\omega\sin\omega\theta
    Split ff into its frequency components.
  2. Step 2:
    fω(θ+s)−fω(θ−s)=2sin⁡ωsω fω′(θ)f_\omega(\theta+s)-f_\omega(\theta-s)=\frac{2\sin\omega s}{\omega}\,f_\omega'(\theta)
    The shift identity of the general rule with 2r=ω2r=\omega, applied to each component separately.
  3. Step 3:
    12[f(θ+π2)−f(θ−π2)]=sin⁡π412f1/2′(θ)+sin⁡π21f1′(θ)=2 f1/2′(θ)+f1′(θ)\tfrac12\Big[f\big(\theta+\tfrac\pi2\big)-f\big(\theta-\tfrac\pi2\big)\Big]=\frac{\sin\frac\pi4}{\frac12}f_{1/2}'(\theta)+\frac{\sin\frac\pi2}{1}f_1'(\theta)=\sqrt2\,f_{1/2}'(\theta)+f_1'(\theta)
    The two-term formula is linear in ff, so it acts on each component with the factor sin⁡(ωπ/2)ω\frac{\sin(\omega\pi/2)}{\omega}.
  4. Step 4:
    f′(θ)=f1/2′(θ)+f1′(θ) ⇒ 12[f(θ+π2)−f(θ−π2)]−f′(θ)=(2−1)f1/2′(θ)f'(\theta)=f_{1/2}'(\theta)+f_1'(\theta)\ \Rightarrow\ \tfrac12\Big[f\big(\theta+\tfrac\pi2\big)-f\big(\theta-\tfrac\pi2\big)\Big]-f'(\theta)=\big(\sqrt2-1\big)f_{1/2}'(\theta)
    The two-term rule overweights the frequency-12\tfrac12 component by the factor 2\sqrt2. It is correct only when that component has zero derivative, for example when the control qubit is in a basis state.
Theorem

Four-term shift rule for controlled rotations

For every controlled Pauli rotation and every θ\theta, f′(θ)=c+[f(θ+π2)−f(θ−π2)]−c−[f(θ+3π2)−f(θ−3π2)],c±=2±142,f'(\theta)=c_+\Big[f\big(\theta+\tfrac\pi2\big)-f\big(\theta-\tfrac\pi2\big)\Big]-c_-\Big[f\big(\theta+\tfrac{3\pi}2\big)-f\big(\theta-\tfrac{3\pi}2\big)\Big],\qquad c_\pm=\frac{\sqrt2\pm1}{4\sqrt2}, i.e. c+=0.4267766953c_+=0.4267766953 and c−=0.0732233047c_-=0.0732233047. The coefficients are the unique solution of a 2×22\times2 linear system expressing exactness on the two frequencies ω∈{12,1}\omega\in\{\tfrac12,1\}.

Proof (8 steps)
  1. Step 1:
    D(θ):=c+[f(θ+α)−f(θ−α)]−c−[f(θ+β)−f(θ−β)],α=π2, β=3π2D(\theta):=c_+\big[f(\theta+\alpha)-f(\theta-\alpha)\big]-c_-\big[f(\theta+\beta)-f(\theta-\beta)\big],\qquad\alpha=\tfrac\pi2,\ \beta=\tfrac{3\pi}2
    Ansatz: a combination of two symmetric differences with unknown coefficients. Symmetric differences remove the constant a0a_0 automatically.
  2. Step 2:
    D(θ)=∑ω∈{12,1}2ω[c+sin⁡ωα−c−sin⁡ωβ] fω′(θ)D(\theta)=\sum_{\omega\in\{\frac12,1\}}\frac{2}{\omega}\big[c_+\sin\omega\alpha-c_-\sin\omega\beta\big]\,f_\omega'(\theta)
    Apply the per-component shift identity fω(θ+s)−fω(θ−s)=2sin⁡ωsωfω′(θ)f_\omega(\theta+s)-f_\omega(\theta-s)=\frac{2\sin\omega s}{\omega}f_\omega'(\theta) to both differences.
  3. Step 3:
    D≡f′  ⟺  2ω[c+sin⁡ωα−c−sin⁡ωβ]=1for ω=12 and ω=1D\equiv f'\iff\frac2\omega\big[c_+\sin\omega\alpha-c_-\sin\omega\beta\big]=1\quad\text{for }\omega=\tfrac12\text{ and }\omega=1
    D=f′D=f' for every controlled-rotation cost function if and only if each frequency component is reproduced with weight 11, because f1/2′f_{1/2}' and f1′f_1' vary independently with the coefficients a1,b1,a2,b2a_1,b_1,a_2,b_2.
  4. Step 4:
    ω=1:  2[c+⋅1−c−⋅(−1)]=1;ω=12:  4[c+22−c−22]=1\omega=1:\ \ 2\big[c_+\cdot1-c_-\cdot(-1)\big]=1;\qquad\omega=\tfrac12:\ \ 4\Big[c_+\tfrac{\sqrt2}2-c_-\tfrac{\sqrt2}2\Big]=1
    sin⁡π2=1\sin\frac\pi2=1, sin⁡3π2=−1\sin\frac{3\pi}2=-1, sin⁡π4=sin⁡3π4=22\sin\frac\pi4=\sin\frac{3\pi}4=\frac{\sqrt2}2.
  5. Step 5:
    (111−1)(c+c−)=(12122),det⁡(111−1)=−2≠0\begin{pmatrix}1&1\\1&-1\end{pmatrix}\begin{pmatrix}c_+\\c_-\end{pmatrix}=\begin{pmatrix}\tfrac12\\\tfrac1{2\sqrt2}\end{pmatrix},\qquad\det\begin{pmatrix}1&1\\1&-1\end{pmatrix}=-2\ne0
    Divide the first equation by 22 and the second by 222\sqrt2. The determinant is non-zero, so the solution exists and is unique.
  6. Step 6:
    c+=12(12+122)=2+142,c−=12(12−122)=2−142c_+=\tfrac12\Big(\tfrac12+\tfrac1{2\sqrt2}\Big)=\frac{\sqrt2+1}{4\sqrt2},\qquad c_-=\tfrac12\Big(\tfrac12-\tfrac1{2\sqrt2}\Big)=\frac{\sqrt2-1}{4\sqrt2}
    Add and subtract the two equations, then put 12±122=2±122\tfrac12\pm\tfrac1{2\sqrt2}=\tfrac{\sqrt2\pm1}{2\sqrt2} on a common denominator.
  7. Step 7:
    c+=1.4142135624+15.6568542495=0.4267766953,c−=0.41421356245.6568542495=0.0732233047c_+=\frac{1.4142135624+1}{5.6568542495}=0.4267766953,\qquad c_-=\frac{0.4142135624}{5.6568542495}=0.0732233047
    Numerical values; they satisfy c++c−=12c_++c_-=\tfrac12 and c+−c−=122=0.3535533906c_+-c_-=\tfrac1{2\sqrt2}=0.3535533906.
  8. Step 8:
    f=a0+f1/2+f1 ⇒ D(θ)=f1/2′(θ)+f1′(θ)=f′(θ)f=a_0+f_{1/2}+f_1\ \Rightarrow\ D(\theta)=f_{1/2}'(\theta)+f_1'(\theta)=f'(\theta)
    By the spectrum theorem every controlled-rotation cost function has only these components, so the rule is exact. It needs the four evaluations f(θ±π2)f(\theta\pm\tfrac\pi2) and f(θ±3π2)f(\theta\pm\tfrac{3\pi}2).

∎

In four-term mode the visualiser draws a controlled-rotation cost function with both frequencies. It marks the four evaluation points θ±π2\theta\pm\tfrac\pi2 and θ±3π2\theta\pm\tfrac{3\pi}2, and compares the four-term estimate with the tangent and with the failing two-term estimate. With the default coefficients (a1=1a_1=1, a2=12a_2=\tfrac12) the frequency-12\tfrac12 component is present, so the two-term value is visibly wrong.

There is an alternative that uses only two-term rules. As a time-ordered circuit, CNOT(c,t)\mathrm{CNOT}(c,t), RY(t)(−θ2)R_Y^{(t)}(-\tfrac\theta2), CNOT(c,t)\mathrm{CNOT}(c,t), RY(t)(θ2)R_Y^{(t)}(\tfrac\theta2) equals CRY(θ)CR_Y(\theta). With the control in ∣0⟩\lvert0\rangle the two rotations cancel. With the control in ∣1⟩\lvert1\rangle, XRY(−θ2)X=RY(θ2)XR_Y(-\tfrac\theta2)X=R_Y(\tfrac\theta2), so the target is rotated by θ2+θ2=θ\tfrac\theta2+\tfrac\theta2=\theta. The parameter θ\theta now occurs twice, with scales a=−12a=-\tfrac12 and a=+12a=+\tfrac12, and the chain rule again costs four evaluations. More generally, Wierichs et al. (2022) show that a cost function with RR distinct positive frequencies admits an exact rule with 2R2R evaluations; the controlled rotation is the case R=2R=2.

Circuit · 2 qubits

Quantum circuit on 2 qubitsX on q1 controlled by q0; RY(-0.40) on q1; X on q1 controlled by q0; RY(0.40) on q1q0q1RY(-0.40)RY(0.40)

CRY(θ)CR_Y(\theta) at θ=0.8\theta=0.8 written with two CNOTs: CNOT(0,1)\mathrm{CNOT}(0,1), RY(−0.4)R_Y(-0.4) on qubit 11, CNOT(0,1)\mathrm{CNOT}(0,1), RY(0.4)R_Y(0.4) on qubit 11. Its unitary equals the controlled rotation ∣0⟩⟨0∣⊗I+∣1⟩⟨1∣⊗RY(0.8)\lvert0\rangle\langle0\rvert\otimes I+\lvert1\rangle\langle1\rvert\otimes R_Y(0.8) exactly.

U=(10000100000.9211−0.3894000.38940.9211)U = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0.9211 & -0.3894 \\ 0 & 0 & 0.3894 & 0.9211 \end{pmatrix}

Cost of a gradient and an implementation#

Proposition

Number of circuit evaluations per gradient

Let a circuit contain mPm_{\mathrm{P}} parameterised Pauli rotations and mCm_{\mathrm{C}} parameterised controlled Pauli rotations. The exact gradient with respect to all parameters by shift rules costs 2mP+4mC2m_{\mathrm{P}}+4m_{\mathrm{C}} circuit evaluations, and the unshifted value f(θ)f(\theta) is not among them. In particular, if each of PP parameters occupies a single Pauli rotation, the cost is exactly 2P2P evaluations, or 2PN2PN shots with NN shots per evaluation. For comparison, forward differences use P+1P+1 evaluations and central differences 2P2P, both with a systematic error.

Proof (6 steps)
  1. Step 1:
    ∂jf=12[f(θ+π2ej)−f(θ−π2ej)] ⇒ 2 evaluations per Pauli occurrence\partial_jf=\tfrac12\Big[f\big(\theta+\tfrac\pi2e_j\big)-f\big(\theta-\tfrac\pi2e_j\big)\Big]\ \Rightarrow\ 2\ \text{evaluations per Pauli occurrence}
    The two-term rule for one occurrence; with repeated parameters the chain rule sums over occurrences, each with its own two evaluations.
  2. Step 2:
    θ+π2ej≠θ±π2ek(j≠k)\theta+\tfrac\pi2e_j\ne\theta\pm\tfrac\pi2e_k\quad(j\ne k)
    The shifted points of different parameters differ in two coordinates, so no evaluation can be shared between parameters.
  3. Step 3:
    {θ±π2eg, θ±3π2eg} ⇒ 4 evaluations per controlled occurrence\Big\{\theta\pm\tfrac\pi2e_g,\ \theta\pm\tfrac{3\pi}2e_g\Big\}\ \Rightarrow\ 4\ \text{evaluations per controlled occurrence}
    The four-term rule for one controlled rotation; the points are distinct because π2\tfrac\pi2 and 3π2\tfrac{3\pi}2 differ modulo 4π4\pi, the period of ff for frequencies 12\tfrac12 and 11.
  4. Step 4:
    #evaluations=2mP+4mC;mP=P, mC=0 ⇒ 2P\#\text{evaluations}=2m_{\mathrm{P}}+4m_{\mathrm{C}};\qquad m_{\mathrm{P}}=P,\ m_{\mathrm{C}}=0\ \Rightarrow\ 2P
    Sum over all occurrences.
  5. Step 5:
    coefficient of f(θ) in every rule=0\text{coefficient of }f(\theta)\text{ in every rule}=0
    All rules are combinations of symmetric differences f(θ+s)−f(θ−s)f(\theta+s)-f(\theta-s), so the unshifted value never enters. A training loop that also needs the loss value spends one extra evaluation for it, which is not part of the gradient.
  6. Step 6:
    forward FD: f(θ+hej)−f(θ)h ⇒ P+1;central FD: 2P;bias O(h) and O(h2)\text{forward FD: }\frac{f(\theta+he_j)-f(\theta)}{h}\ \Rightarrow\ P+1;\qquad\text{central FD: }2P;\qquad\text{bias }O(h)\text{ and }O(h^2)
    Forward differences share f(θ)f(\theta) between all parameters; both difference schemes are approximations with a step-size bias, unlike the shift rules.

∎

Quick check

An ansatz has 66 parameters in single Pauli rotations and 22 further parameters, each in a single controlled rotation CRYCR_Y. How many circuit evaluations does one exact gradient by shift rules cost?

tolerance ±0.5

The TypeScript cell implements the rule with the simulator's exact expectationValue on a three-qubit circuit with five parameters. It compares the result with central finite differences, counts 2P=102P=10 evaluations, and shows that the general rule with r=12r=\tfrac12 returns the same derivative for every shift ss. Finally, on a CRYCR_Y gate it compares the failing two-term rule, the four-term rule and finite differences.

TypeScriptrunnable
import { FOUR_TERM_C_MINUS, FOUR_TERM_C_PLUS, expectationValue, finiteDifferenceGradient, mulberry32 } from '@/lib/quantum-simulator';

// A three-qubit circuit in which every parameter sits in one Pauli rotation R_P(theta) = exp(-i theta P / 2).
const names = ['a', 'b', 'c', 'd', 'e'];
const rot = (gate: string, q: number, name: string) => ({ kind: 'gate', gate, targets: [q], params: [{ param: name }] });
const cnot = (c: number, t: number) => ({ kind: 'gate', gate: 'CNOT', targets: [c, t] });
const circuit: any = {
  version: 1, numQubits: 3, numClbits: 0,
  parameters: names.map((name) => ({ name, value: 0 })),
  ops: [rot('RY', 0, 'a'), rot('RX', 1, 'b'), cnot(0, 1), rot('RY', 1, 'c'), cnot(1, 2), rot('RX', 2, 'd'), cnot(2, 0), rot('RY', 0, 'e')],
};
const H = [{ coeff: 1, ops: 'ZIZ' }, { coeff: -0.7, ops: 'XYX' }, { coeff: 0.5, ops: 'IZI' }];

let calls = 0;
const f = (theta: number[]): number => {
  calls++;
  const bindings: Record<string, number> = {};
  names.forEach((n, i) => (bindings[n] = theta[i]));
  return expectationValue(circuit, H, bindings);
};
const shift = (theta: number[], j: number, s: number) => theta.map((v, k) => (k === j ? v + s : v));

// Parameter-shift gradient: two evaluations per parameter, shifts +-pi/2, prefactor 1/2.
function psrGradient(theta: number[]): number[] {
  return theta.map((_, j) => 0.5 * (f(shift(theta, j, Math.PI / 2)) - f(shift(theta, j, -Math.PI / 2))));
}

const rng = mulberry32(48);
const theta = names.map(() => 2 * Math.PI * rng());
const g = psrGradient(theta);
const bindings = Object.fromEntries(names.map((n, i) => [n, theta[i]]));
const fd = finiteDifferenceGradient(circuit, H, { bindings }).gradient;
console.log('theta      : ' + theta.map((v) => v.toFixed(4)).join(' '));
console.log('PSR        : ' + g.map((v) => v.toFixed(8)).join(' '));
console.log('central FD : ' + names.map((n) => fd[n].toFixed(8)).join(' '));
const dev = Math.max(...names.map((n, i) => Math.abs(g[i] - fd[n])));
console.log('max |PSR - FD| < 1e-8: ' + (dev < 1e-8) + ', evaluate calls: ' + calls + ' = 2P');

// The general rule r / sin(2 r s) [f(theta + s) - f(theta - s)] with r = 1/2 gives the same value for every shift s.
for (const s of [0.1, 1.0, Math.PI / 2, 2.5]) {
  const d = (0.5 / Math.sin(s)) * (f(shift(theta, 0, s)) - f(shift(theta, 0, -s)));
  console.log('s = ' + s.toFixed(4) + ': df/da = ' + d.toFixed(8));
}

// Controlled rotation CRY(t): generator eigenvalues {0, 0, 1/2, -1/2}, frequencies 1/2 and 1.
const ctrl: any = {
  version: 1, numQubits: 2, numClbits: 0, parameters: [{ name: 't', value: 0 }],
  ops: [
    { kind: 'gate', gate: 'RY', targets: [0], params: [1.1] },
    { kind: 'gate', gate: 'RX', targets: [1], params: [0.4] },
    { kind: 'gate', gate: 'CRY', targets: [0, 1], params: [{ param: 't' }] },
  ],
};
const H2 = [{ coeff: 1, ops: 'XZ' }, { coeff: 0.5, ops: 'IZ' }];
const g2 = (t: number) => expectationValue(ctrl, H2, { t });
const t0 = 0.8;
const twoTerm = 0.5 * (g2(t0 + Math.PI / 2) - g2(t0 - Math.PI / 2));
const fourTerm = FOUR_TERM_C_PLUS * (g2(t0 + Math.PI / 2) - g2(t0 - Math.PI / 2))
  - FOUR_TERM_C_MINUS * (g2(t0 + 1.5 * Math.PI) - g2(t0 - 1.5 * Math.PI));
const fdT = finiteDifferenceGradient(ctrl, H2, { bindings: { t: t0 } }).gradient.t;
console.log('CRY: two-term ' + twoTerm.toFixed(8) + ', four-term ' + fourTerm.toFixed(8) + ', central FD ' + fdT.toFixed(8));

Expected output

theta      : 5.3106 3.3600 1.5749 0.4384 2.9134
PSR        : -0.58333288 0.49028787 0.27549995 0.06730340 0.51019828
central FD : -0.58333288 0.49028787 0.27549995 0.06730340 0.51019828
max |PSR - FD| < 1e-8: true, evaluate calls: 10 = 2P
s = 0.1000: df/da = -0.58333288
s = 1.0000: df/da = -0.58333288
s = 1.5708: df/da = -0.58333288
s = 2.5000: df/da = -0.58333288
CRY: two-term -0.31628747, four-term -0.25008444, central FD -0.25008444

The NumPy cell repeats the checks with explicit 4×44\times4 matrices: the two-term rule against finite differences, the coefficients c±c_\pm from the 2×22\times2 system, the four-term rule on CRXCR_X, CRYCR_Y and CRZCR_Z, and the repeated-parameter example, where the naive global shift returns 00.

Pythonrunnable
import numpy as np

I2 = np.eye(2)
X = np.array([[0, 1], [1, 0]], dtype=complex)
Y = np.array([[0, -1j], [1j, 0]])
Z = np.diag([1.0 + 0j, -1.0])
P0, P1 = np.diag([1.0, 0.0]), np.diag([0.0, 1.0])


def rot(P, t):  # R_P(t) = exp(-i t P / 2) = cos(t/2) I - i sin(t/2) P
    return np.cos(t / 2) * I2 - 1j * np.sin(t / 2) * P


def crot(P, t):  # control qubit 0, target qubit 1 (qubit 0 = most significant)
    return np.kron(P0, I2) + np.kron(P1, rot(P, t))


CNOT = np.kron(P0, I2) + np.kron(P1, X)
psi0 = np.array([1, 0, 0, 0], dtype=complex)

# (1) Two-term rule on R_Y(a) (x) R_X(b), CNOT, observable H = Z1 + 0.5 X0 X1.
H = np.kron(I2, Z) + 0.5 * np.kron(X, X)


def f(a, b):
    psi = CNOT @ np.kron(rot(Y, a), rot(X, b)) @ psi0
    return np.real(np.vdot(psi, H @ psi))


a, b, h = 0.9, -0.4, 1e-6
psr = [0.5 * (f(a + np.pi / 2, b) - f(a - np.pi / 2, b)), 0.5 * (f(a, b + np.pi / 2) - f(a, b - np.pi / 2))]
fd = [(f(a + h, b) - f(a - h, b)) / (2 * h), (f(a, b + h) - f(a, b - h)) / (2 * h)]
print("two-term PSR :", " ".join(f"{g:+.8f}" for g in psr))
print("central FD   :", " ".join(f"{g:+.8f}" for g in fd))
print("agree to 1e-8:", bool(np.max(np.abs(np.subtract(psr, fd))) < 1e-8))

# (2) Controlled rotations: frequencies 1/2 and 1, so the two-term rule fails and the four-term rule is exact.
cp = (np.sqrt(2) + 1) / (4 * np.sqrt(2))
cm = (np.sqrt(2) - 1) / (4 * np.sqrt(2))
coef = np.linalg.solve(np.array([[1.0, 1.0], [1.0, -1.0]]), np.array([0.5, 1 / (2 * np.sqrt(2))]))
print(f"c+ = {cp:.10f}  c- = {cm:.10f}  (2x2 solve: {coef[0]:.10f}, {coef[1]:.10f})")
H2 = np.kron(X, Z) + 0.5 * np.kron(I2, Z)
prep = np.kron(rot(Y, 1.1), rot(X, 0.4))
for name, P in (("CRX", X), ("CRY", Y), ("CRZ", Z)):
    def g(t):
        psi = crot(P, t) @ np.kron(I2, rot(Y, 0.7)) @ prep @ psi0
        return np.real(np.vdot(psi, H2 @ psi))
    t = 0.8
    two = 0.5 * (g(t + np.pi / 2) - g(t - np.pi / 2))
    four = cp * (g(t + np.pi / 2) - g(t - np.pi / 2)) - cm * (g(t + 3 * np.pi / 2) - g(t - 3 * np.pi / 2))
    exact = (g(t + h) - g(t - h)) / (2 * h)
    print(f"{name}: two-term {two:+.8f}  four-term {four:+.8f}  central FD {exact:+.8f}")

# (3) A parameter used twice: R_Y(alpha) R_Y(alpha) |0> = R_Y(2 alpha) |0>, F = <Z> = cos(2 alpha).
F = lambda al: np.cos(2 * al)
two_occ = lambda t1, t2: np.cos(t1 + t2)  # the same circuit with independent occurrence angles
al = 0.6
naive = 0.5 * (F(al + np.pi / 2) - F(al - np.pi / 2))
chain = sum(0.5 * (two_occ(*(np.array([al, al]) + np.pi / 2 * e)) - two_occ(*(np.array([al, al]) - np.pi / 2 * e)))
            for e in np.eye(2))
print(f"repeated parameter: naive shift {round(naive, 12) + 0.0:+.8f}  sum over occurrences {chain:+.8f}  exact {-2 * np.sin(2 * al):+.8f}")

Expected output

two-term PSR : -0.41068688 +0.24206632
central FD   : -0.41068688 +0.24206632
agree to 1e-8: True
c+ = 0.4267766953  c- = 0.0732233047  (2x2 solve: 0.4267766953, 0.0732233047)
CRX: two-term -0.50500255  four-term -0.38816465  central FD -0.38816465
CRY: two-term -0.64278888  four-term -0.49127925  central FD -0.49127925
CRZ: two-term -0.17287827  four-term -0.12224339  central FD -0.12224339
repeated parameter: naive shift +0.00000000  sum over occurrences -1.86407817  exact -1.86407817

The lab asks for exactly this function. It receives a black-box evaluate and must return the gradient using precisely 2P2P calls; the grader counts the calls.

References

  1. K. Mitarai, M. Negoro, M. Kitagawa and K. Fujii, Quantum circuit learning, Physical Review A 98, 032309 (2018).doi:10.1103/PhysRevA.98.032309
  2. 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
  3. D. Wierichs, J. Izaac, C. Wang and C. Y.-Y. Lin, General parameter-shift rules for quantum gradients, Quantum 6, 677 (2022).doi:10.22331/q-2022-03-30-677
  4. G. E. Crooks, Gradients of parameterized quantum gates using the parameter-shift rule and gate decomposition, arXiv:1905.13311 (2019).arxiv.org
  5. G.-L. R. Anselmetti, D. Wierichs, C. Gogolin and R. M. Parrish, Local, expressive, quantum-number-preserving VQE ansätze for fermionic systems, New Journal of Physics 23, 113010 (2021).doi:10.1088/1367-2630/ac2cb3

Hands-on labs

  • Gradients with the parameter-shift rule

    1 step · auto-graded

    Implement the parameter-shift gradient of a variational circuit using only a black-box evaluate(params) that returns the exact expectation value f(θ)=⟨ψ(θ)∣H∣ψ(θ)⟩f(\theta)=\langle\psi(\theta)\vert H\vert\psi(\theta)\rangle. Every parameter sits in a single Pauli rotation RP(θj)=e−iθjP/2R_P(\theta_j)=e^{-i\theta_jP/2}, so ∂jf(θ)=12[f(θ+π2ej)−f(θ−π2ej)]\partial_jf(\theta)=\tfrac12\big[f(\theta+\tfrac\pi2e_j)-f(\theta-\tfrac\pi2e_j)\big] is exact. The grader checks the values on one visible and four hidden random three-qubit circuits, separately checks that your function calls evaluate exactly 2P2P times, and finally checks that the values are exact to 10−1210^{-12}, which no finite-difference quotient achieves.

    Open lab