Chapter 16: PennyLane

"Qiskit gives you a measurement. PennyLane gives you a derivative. Almost everything else follows from that choice."

Overview

Q# added ceremony to catch errors. PennyLane removes it to make circuits differentiable.

The organizing idea is simple and consequential: a quantum circuit is a function from parameters to an expectation value, and if you can differentiate it, you can drop it into any machine-learning optimizer. PennyLane makes the gradient the primary output rather than the measurement.

And it does so exactly. Not by finite differences — by the parameter-shift rule, which computes an exact derivative from two evaluations of the same circuit at shifted parameter values. It works on hardware, where you cannot take an infinitesimal step, and it is exact rather than approximate:

$$\frac{\partial}{\partial\theta}\langle Z\rangle = \frac{f(\theta + \pi/2) - f(\theta - \pi/2)}{2}$$

Verified against the analytic answer to machine precision in §16.3.

That machinery makes variational quantum algorithms practical, and it leads directly to the reason they are hard. Measure how large a typical gradient is as you add qubits:

   qubits    Var[gradient]    max|gradient|
        2       1.0431e-01           0.8841
        6       8.0048e-03           0.3712
       10       4.6517e-04           0.0920

   fit: Var ~ exp(-0.676 n)     ->  x0.51 per qubit added

The gradient variance halves with every qubit you add. Over eight qubits it fell by 224×. Extrapolated to fifty qubits, resolving a gradient would take roughly 10¹⁵ shots.

This is the barren plateau, and it is the central obstacle to scaling variational quantum algorithms. Chapter 24 will build VQE properly; this chapter measures the wall it runs into.

In this chapter, you will learn to:

• Build a QNode and differentiate it.
• Derive and verify the parameter-shift rule.
• Count what a gradient costs.
• Run and diagnose an optimization loop.
• Measure the barren plateau.
• Tell a flat landscape from a structurally zero gradient.

Learning Paths

How to read this chapter by track. - 🔰 Beginner — §16.2 and §16.3. The parameter-shift rule is the idea to take away. - 🔬 Researcher — §16.6 and §16.7. Barren plateaus decide what is worth attempting. - 🤖 Quantum ML — all of it; this is your chapter and your framework. - 🏗️ Quantum Engineer — §16.4's cost accounting, which is what your shot budget is spent on. - 🔐 Security — skim; variational methods are largely orthogonal to your track.

16.1 The QNode

import pennylane as qml
from pennylane import numpy as pnp

dev = qml.device("default.qubit", wires=1)

@qml.qnode(dev)
def circuit(theta):
    qml.RY(theta, wires=0)
    return qml.expval(qml.PauliZ(0))

A QNode is a quantum function bound to a device. The decorator turns an ordinary Python function into something that runs a circuit and — critically — can be differentiated.

Three things differ from everything so far.

There is no circuit object. You do not build a QuantumCircuit and pass it somewhere; the function body is the circuit, executed each time you call it. This is why parameters are plain Python arguments rather than Parameter objects (Qiskit) or sympy symbols (Cirq) — and why Chapter 8's parameter-ordering bug cannot occur here either, for a third distinct reason.

The return statement declares the measurement. qml.expval, qml.probs, qml.sample, qml.state — what you return determines what the device computes. Returning an expectation value is the default because expectation values are what you differentiate.

pennylane.numpy wraps NumPy to track gradients. Arrays carry requires_grad, exactly as in an autodiff framework, because that is what it is.

  <Z>(0.7)  = 0.764842        analytic cos(0.7) = 0.764842

wires rather than qubits is deliberate: PennyLane also supports continuous-variable devices, where "qubit" would be wrong.

What the abstraction costs

Nothing in this book is free, and the QNode's convenience is bought with something specific: there is no artifact.

In Qiskit a QuantumCircuit is an object. You can hold it, print it, transpile it (Chapter 10), serialize it to QASM and hand it to another framework (Chapter 6), count its gates before you run it, and pass it to three different backends. It exists independently of any execution.

A QNode has none of that, because the circuit is a side effect of calling the function. Every inspection tool therefore has to run the trace to get anything at all: qml.draw(qnode)(args) calls the function, qml.specs(qnode)(args) calls the function, and qml.workflow.construct_tape(qnode)(args) calls the function to hand you the tape it built. You cannot look at the circuit without supplying arguments, because until you do there is no circuit — only Python that would produce one.

There is a second, subtler cost. qml.specs reports what the QNode looks like at the level you ask for, and the default is not the gate list:

   Total gates: 1
   Gate counts:
     - StronglyEntanglingLayers: 1
   Depth: 1

A six-layer, 54-parameter ansatz reported as one gate of depth one. That is not wrong — the template is one operation until something decomposes it — but a practitioner reading "depth 1" and concluding the circuit is shallow has read the abstraction rather than the circuit. Chapter 10 §10.7's insistence on reading the diff — the gate counts before and after each stage, rather than whichever count came to hand — applies here with a PennyLane accent.

The trade is Chapter 14 §14.1's, seen from the other end. Cirq gives you less abstraction and fewer surprises; PennyLane gives you more abstraction and one specific surprise — that the object you are reasoning about does not exist between calls.

⚠️ Common Pitfall: plain NumPy produces no error and no gradient.

This is the single most common PennyLane bug in this book, it costs about twenty minutes the first time, and it is two characters wide. Chapter 32 lost time to it; the instructor guide lists it as struggle number 8.

qml.grad differentiates with respect to trainable arguments, and trainability is a property of the array, not of the QNode. Measured on §16.1's circuit at $\theta = 0.7$:

```text argument circuit value gradient shape pnp.array(0.7, requires_grad=True) 0.764842 -0.644218 () pnp.array(0.7) [pnp default] 0.764842 -0.644218 () pnp.array(0.7, requires_grad=False) 0.764842 () (0,) np.array(0.7) [plain numpy] 0.764842 () (0,) 0.7 [Python float] 0.764842 () (0,)

plain numpy vector [0.7, 0.3] -> array([], dtype=float64) shape (0,) pnp vector, requires_grad=True -> [-0.61544466, -0.22602632] shape (2,) ```

Read the middle column. The circuit value is 0.764842 in every row — correct, identical, and completely uninformative about whether the gradient will work. The failure is invisible in the forward pass and only appears in the derivative.

And it is a shape, not an exception. You get back an empty tuple or an empty array, which then flows into your optimizer, which updates nothing, which produces a cost curve that is perfectly flat. A flat cost curve is exactly what §16.6's barren plateau also produces — so the cheapest bug in the chapter and the deepest obstacle in the field present with the same symptom. Check the gradient's shape before you conclude anything about a landscape.

Note row two: pnp.array(0.7) works. pennylane.numpy defaults requires_grad to True, so the keyword is documentation rather than magic. The bug is importing the wrong NumPy, and it is almost always a stray import numpy as np at the top of a file that also does from pennylane import numpy as pnp, followed by one array built with the wrong one.

PennyLane 0.45.1 does raise a UserWarning — "Attempted to differentiate a function with no trainable parameters" — which is genuinely helpful and genuinely easy to miss in a loop that prints a cost every ten steps. Promote it: warnings.simplefilter("error", UserWarning) during development turns twenty minutes into one traceback.

🗝️ Version Note: PennyLane 0.45.1, and where shots now live.

Two API details in this chapter have moved, and both bite silently.

Device-level shots are deprecated. The idiom every tutorial written before 2025 uses —

python dev = qml.device("default.qubit", wires=1, shots=1000)

— still runs in 0.45.1 but emits PennyLaneDeprecationWarning: Setting shots on device is deprecated. Please use the set_shots transform on the respective QNode instead. The current form puts shots on the QNode, where they belong, because shots are a property of how you execute rather than of the hardware:

python @qml.set_shots(1000) @qml.qnode(dev) def circuit(theta): ...

This matters more than a deprecation usually does, because shots determine which differentiation methods are legal — see §16.3's diff_method table. Moving shots from the device to the QNode moves that decision to the same place.

The workflow introspection namespace. qml.workflow.construct_tape(qnode)(args) returns the tape a call would build, and qml.workflow.get_best_diff_method(qnode)(args) returns the method diff_method="best" will actually resolve to. Neither is in the top-level qml namespace. Both are the fastest way to answer "what is this QNode really doing," and §16.3 uses the second one to establish something the chapter would otherwise have got wrong.

diff_method=None is not "the default." It means not differentiable and raises QuantumFunctionError: Derivatives cannot be calculated with diff_method=None the moment you call qml.grad. The default is "best".

16.2 Differentiation

theta = pnp.array(0.7, requires_grad=True)
gradient = qml.grad(circuit)(theta)

  d<Z>/dtheta = -0.644218      analytic -sin(0.7) = -0.644218
  error       = 0.00e+00

Exact to machine precision. For $\langle Z\rangle = \cos\theta$ the derivative is $-\sin\theta$, and PennyLane returns it with zero error.

That is worth pausing on, because the obvious implementation — finite differences — cannot do this.

16.3 The Parameter-Shift Rule

The problem. On hardware you cannot evaluate a derivative by taking an infinitesimal step. Each evaluation costs thousands of shots and carries statistical noise, so a finite difference divides a small, noisy number by a small number.

The solution. For a gate of the form $e^{-i\theta P/2}$ with $P$ a Pauli, the expectation value is exactly sinusoidal in $\theta$. A sinusoid's derivative can be recovered from two samples a quarter period apart:

$$\frac{\partial f}{\partial\theta} = \frac{f(\theta + \tfrac{\pi}{2}) - f(\theta - \tfrac{\pi}{2})}{2}$$

Not an approximation. Verified by hand:

s = np.pi / 2
manual = (circuit(0.7 + s) - circuit(0.7 - s)) / 2

  manual parameter shift = -0.644218
  qml.grad               = -0.644218
  analytic               = -0.644218

Compare with finite differences on the same circuit:

  h=1e-01:  -0.643144528    error 1.07e-03
  h=1e-03:  -0.644217580    error 1.07e-07
  h=1e-05:  -0.644217687    error 1.98e-11
  h=1e-07:  -0.644217688    error 2.91e-10

Finite differences have a sweet spot and get worse on both sides of it. Too large and truncation error dominates; too small and floating-point cancellation does — visible above as the error rising again from $10^{-5}$ to $10^{-7}$. There is no step size that gives an exact answer, and on hardware the noise floor moves the sweet spot to a much worse place.

The parameter-shift rule has no step size to tune. $\pi/2$ is not small; it is a large shift, which is precisely why the rule survives on noisy hardware where finite differences do not.

📐 Math Aside — why $\pi/2$, and why this only works for some gates.

For a single parameterized rotation $U(\theta) = e^{-i\theta P/2}$ with $P^2 = I$, any expectation value takes the form

$$f(\theta) = A\cos\theta + B\sin\theta$$

Then $f'(\theta) = -A\sin\theta + B\cos\theta$, and evaluating at $\theta \pm \pi/2$:

$$f(\theta + \tfrac{\pi}{2}) - f(\theta - \tfrac{\pi}{2}) = 2\left(B\cos\theta - A\sin\theta\right) > = 2f'(\theta)$$

The rule is exact because the function is exactly a sinusoid — a consequence of $P$ having eigenvalues $\pm 1$.

It does not apply to every gate. A gate whose generator has more than two distinct eigenvalues produces a function with more frequencies, requiring more shift points (a four-term rule, or a general decomposition). PennyLane handles this per gate; diff_method="parameter-shift" will tell you if it cannot.

📐 Math Aside — where $A$ and $B$ actually come from, and why $\pi/2$ rather than $0.1$.

The aside above asserts the sinusoidal form. Deriving it takes four lines and pays for itself, because the derivation says exactly which gates are covered and exactly why the shift is what it is.

Because $P^2 = I$, the exponential closes in two terms:

$$U(\theta) = e^{-i\theta P/2} = \cos\tfrac{\theta}{2}\,I - i\sin\tfrac{\theta}{2}\,P$$

Let $|\psi\rangle$ be the state entering the gate and $O$ the observable propagated back to just after it. Then $f(\theta) = \langle\psi|U^\dagger O U|\psi\rangle$, and expanding:

$$U^\dagger O U = \cos^2\tfrac{\theta}{2}\,O + \sin^2\tfrac{\theta}{2}\,POP > + \tfrac{i}{2}\sin\theta\,[P, O]$$

Substituting $\cos^2(\theta/2) = (1+\cos\theta)/2$ and $\sin^2(\theta/2) = (1-\cos\theta)/2$ and collecting terms:

$$f(\theta) = C + A\cos\theta + B\sin\theta, \qquad > C = \left\langle\tfrac{O + POP}{2}\right\rangle,\; > A = \left\langle\tfrac{O - POP}{2}\right\rangle,\; > B = \tfrac{i}{2}\big\langle[P,O]\big\rangle$$

Note the constant $C$, which the compact form above suppresses. It is not always zero — it happens to vanish for §16.1's circuit, where $\langle Z\rangle = \cos\theta$ exactly. Fit a circuit where it does not:

text H(0); CRY(theta, [0,1]); measure <Z_1> fitted f(theta) = +0.500000 +0.500000 cos(theta) -0.000000 sin(theta) max residual over 13 points: 4.44e-16

The rule survives it, and that is the point of taking a difference. $C$ appears identically at $\theta + s$ and $\theta - s$, so the subtraction removes it — measured on that circuit at $\theta = 0.7$: the two-term rule returns $-0.322108844$, and $-\tfrac12\sin(0.7) = -0.322108844$. A rule built on a ratio rather than a difference would not have this property.

Now the shift. For any $s$, not just $\pi/2$:

$$f(\theta + s) - f(\theta - s) = 2\sin(s)\big(B\cos\theta - A\sin\theta\big) = 2\sin(s)\,f'(\theta)$$

$$\boxed{\;f'(\theta) = \frac{f(\theta + s) - f(\theta - s)}{2\sin s}\;}$$

Every shift is exact. Measured on §16.1's circuit at $\theta = 0.7$, analytic $-0.644217687$:

text shift s value estimate error 1/(2 sin s) pi/2 1.570796 -0.644217687 2.22e-16 0.5000 pi/3 1.047198 -0.644217687 1.11e-16 0.5774 pi/4 0.785398 -0.644217687 1.11e-16 0.7071 pi/6 0.523599 -0.644217687 2.22e-16 1.0000 0.1 0.100000 -0.644217687 7.77e-16 5.0083 3 pi/4 2.356194 -0.644217687 0.00e+00 0.7071

A shift of $0.1$ is exactly as correct as $\pi/2$. So $\pi/2$ is not chosen for accuracy — it is chosen for noise. Look at the last column: the prefactor $1/(2\sin s)$ multiplies whatever statistical error the two evaluations carry, and it is minimised at $s = \pi/2$, where it equals $1/2$. At $s = 0.1$ it is $5.0083$ — ten times the noise amplification for an identical exact answer.

That is the whole story in one line. A small shift is a finite difference that happens to be exact; $\pi/2$ is the shift that also happens to be optimal. §16.3's noise report below measures what that factor is worth on hardware.

Where the finite-difference error actually comes from

The finite-difference table above is usually read as "small $h$ good, tiny $h$ bad." The real structure is sharper than that, and it can be derived exactly from numbers already on this page.

A central difference on $f(\theta) = C + A\cos\theta + B\sin\theta$ gives

$$\frac{f(\theta + h) - f(\theta - h)}{2h} = \frac{2\sin(h)\,f'(\theta)}{2h} = f'(\theta)\cdot\frac{\sin h}{h}$$

The truncation error is not "some $\mathcal{O}(h^2)$ term." It is exactly the factor $\sin(h)/h$. Which means the entire measured table is predictable to nine digits without running anything:

      h     sin(h)/h          predicted        measured (16.3)     predicted err   measured err
   1e-01   0.998334166     -0.643144528     -0.643144528              1.07e-03       1.07e-03
   1e-03   0.999999833     -0.644217580     -0.644217580              1.07e-07       1.07e-07
   1e-05   0.999999999983  -0.644217687     -0.644217687              1.07e-11       1.98e-11

Two of the three rows match digit for digit, which is a satisfying check on both the derivation and the measurement. Expanding $\sin h/h \approx 1 - h^2/6$ gives the familiar quadratic law, and predicts a truncation error of $|f'|h^2/6$ — which is $1.074\times10^{-3}$ at $h = 10^{-1}$ and $1.074\times10^{-7}$ at $10^{-3}$, both correct.

And the third row is where the story is. Truncation predicts $1.07\times10^{-11}$; the measurement found $1.98\times10^{-11}$. The measurement is worse than truncation theory allows, by roughly a factor of two, because at $h = 10^{-5}$ floating-point cancellation has become the larger of the two errors. That row is the crossover, located precisely: truncation falls as $h^2$, roundoff rises as $1/h$, and the sweet spot sits where the two curves cross. Everything the chapter said qualitatively about "a sweet spot with degradation on both sides" is this crossover, and it is now a number rather than a shape.

⚛️ The Physics Underneath: the shift rule is a statement about the generator's spectrum.

The reason a two-term rule works for RY and not for CRY has nothing to do with PennyLane. It is a fact about eigenvalue gaps.

A gate $e^{-i\theta G}$ makes any expectation value a trigonometric polynomial whose frequencies are the differences between eigenvalues of $G$. For a Pauli generator $P/2$ the eigenvalues are $\pm\tfrac12$, there is exactly one nonzero gap, and $f$ has exactly one frequency. One frequency, two unknowns beyond the constant, two evaluations.

CRY is $|1\rangle\langle1|\otimes R_Y(\theta)$, so its generator is $|1\rangle\langle1|\otimes Y/2$ with eigenvalues $\{0,\,0,\,+\tfrac12,\,-\tfrac12\}$. The distinct gaps are $\tfrac12$ and $1$ — two frequencies, four unknowns, four evaluations. PennyLane knows this per operation, and the execution counts show it doing the accounting:

text one parameter, diff_method="parameter-shift" RY -> 3 executions (2 shifts + 1 forward pass) CRY -> 5 executions (4 shifts + 1 forward pass)

This is why the cost model in §16.4 is stated as $2n+1$ for Pauli rotations and why a hardware ansatz built from controlled rotations is more expensive to differentiate than its gate count suggests. It is also Chapter 22's phase-estimation intuition wearing different clothes: eigenvalue gaps determine how many samples you need to reconstruct a periodic function.

When two terms are not enough

The failure is worth seeing, because it is silent. Take a circuit with a CRY sandwiched between Hadamards on the control, so the half-frequency survives into the measurement:

qml.Hadamard(0)
qml.CRY(theta, wires=[0, 1])
qml.Hadamard(0)
return qml.expval(qml.PauliZ(0))

Sample it, and the period gives the game away:

   f(0)    =  1.000000
   f(pi)   =  0.000000
   f(2 pi) = -1.000000        <- not back where it started
   f(4 pi) =  1.000000        <- period is 4 pi, not 2 pi

A function with period $4\pi$ has a half-frequency component, and no combination of two samples a quarter-period apart can see it. Fitting confirms which basis is right:

   fit to  c + A cos(t) + B sin(t)                      max residual 1.111e+00
   fit to  c + A cos(t/2) + B sin(t/2) + C cos t + D sin t   max residual 4.441e-16

Now differentiate it three ways at $\theta = 0.7$:

   reference (exact simulator)        -0.171448904
   PennyLane parameter-shift          -0.171448904     5 executions
   naive two-term rule                -0.242465365     error 7.10e-02

The naive rule is 41% high, and it does not look wrong. It returns a smooth, plausible, reproducible number of the right sign and order of magnitude — which would produce an optimizer that descends confidently along a systematically incorrect direction. Nothing crashes; the run just converges somewhere else.

The lesson generalises past this gate. If you hand-roll a shift rule — and people do, to save the forward pass or to fuse it into a custom loop — you have quietly assumed something about every generator in your circuit. qml.gradients.param_shift already knows the right rule per operation. This is the book's recurring warning in a new setting: the number that is easy to get is not the number that answers the question, and here the easy number is off by 41% while looking fine.

Which differentiation method did you actually get?

Here is something §16.2 quietly did not tell you, and it changes how to read that result.

qml.grad does not mean "parameter-shift." The QNode's default is diff_method="best", and "best" is resolved against the device:

   device                         diff_method="best" resolves to
   default.qubit, shots=None      backprop
   default.qubit, shots=1000      parameter-shift

§16.2's machine-precision gradient came from backpropagation through the simulator's state vector, not from the parameter-shift rule. §16.3's manual verification is what establishes the rule works; the qml.grad call in §16.2 agreed with it because both are exact, not because they are the same computation.

The costs are not remotely comparable. One gradient of a 4-parameter circuit:

   diff_method        executions   agrees with parameter-shift
   backprop                    1   —
   adjoint                     1   1.67e-16
   finite-diff                 6   5.05e-09
   parameter-shift             9   0.00e+00
   best (no shots)             1   1.67e-16

And the scaling, which is the part that matters:

      n   parameter-shift   finite-diff   backprop      2n+1    n+2
      1                 3             3          1         3      3
      2                 5             4          1         5      4
      4                 9             6          1         9      6
      8                17            10          1        17     10
     16                33            18          1        33     18

Backprop is $\mathcal{O}(1)$ in circuit executions at any parameter count, which is why simulator experiments in Part VI run at all. Finite differences cost $n+2$ under PennyLane's default forward-difference settings — cheaper than parameter-shift's $2n+1$, which is the trap: the method that is wrong on hardware is also the one that looks cheaper on paper.

And then the constraint that decides everything:

   with shots=1000 on default.qubit:
     backprop         -> QuantumFunctionError: device does not support backprop
     adjoint          -> QuantumFunctionError: device does not support adjoint
     parameter-shift  -> 9 executions

Backprop and adjoint differentiation are simulator-only, and not by omission. They require access to the state vector — backprop differentiates through the simulation's arithmetic, adjoint runs the circuit backwards from the final state. A quantum computer has neither. You cannot read out a state, and you cannot run a measurement backwards. The moment shots is set, PennyLane is modelling a device that samples, and both methods become unavailable even on a simulator, which is exactly the right behaviour: it makes the simulator refuse to do something the hardware could not.

📊 What the Numbers Say: a machine-precision gradient on a simulator is not evidence that the parameter-shift rule works.

This is a small thing that misleads reliably. A student runs §16.2, sees error = 0.00e+00, and concludes the parameter-shift rule is exact. The conclusion is true. The evidence does not support it, because the number came from backprop.

The chapter's own structure is the fix: §16.3 computes the shift by hand, with np.pi/2 written explicitly, and compares three ways. That comparison is the evidence; the qml.grad call is not.

It generalises to a habit worth having. When a framework picks a strategy for you, the result validates the framework's choice, not the strategy you had in mind. Chapter 34 §34.4 found the same shape in a library default — zz_feature_map's shipped reps=2 is the worst row of that chapter's own sweep, at 0.6167 against a tuned 0.8500 — and Chapter 29 §29.3 measured a hardware-aware level-1 transpilation at 0.9116 beating a naive level-3 one at 0.7720. A default is a default, not a recommendation.

If you want to know whether parameter-shift works on your ansatz, ask for it by name: @qml.qnode(dev, diff_method="parameter-shift"). If you want to know what it costs, wrap it in a qml.Tracker as §16.4 does. Do not let "best" answer a question you asked about a specific method.

Why exactness matters under shot noise

The chapter has claimed twice now that the parameter-shift rule "survives on noisy hardware where finite differences do not." That is the sort of claim this book is supposed to measure, so measure it.

First derive what to expect. An expectation value of a Pauli from $N$ shots has standard error $\sigma = \sqrt{1 - \langle O\rangle^2}/\sqrt{N}$. Both estimators are a difference of two independent evaluations, so both inherit $\sqrt{2}\,\sigma$ in the numerator — and then they divide by different things:

$$\text{std}\big[\hat{f}'_{\text{shift}}\big] = \frac{\sqrt2\,\sigma}{2}, \qquad \text{std}\big[\hat{f}'_{\text{FD}}\big] = \frac{\sqrt2\,\sigma}{2h}$$

$$\Rightarrow\quad \frac{\text{std}[\hat{f}'_{\text{FD}}]}{\text{std}[\hat{f}'_{\text{shift}}]} \approx \frac{1}{h}$$

The ratio is $1/h$. Not a constant, not a modest penalty — the finite difference is worse by the reciprocal of its own step size, which is the number it was told to make small. The step size that controls truncation error is the same step size that amplifies shot noise, and it moves them in opposite directions. There is no value of $h$ that is good at both, which is a stronger statement than "there is a sweet spot": under shot noise the sweet spot is not merely worse, it has moved to a different place entirely.

📉 Noise Report: the same two gradients at 1,000 shots.

§16.1's circuit at $\theta = 0.7$, analytic derivative $-0.644217687$. Every evaluation drawn from 1,000 shots; 200 independent repeats of each estimator; default.qubit with sampling on.

text method mean std RMS error vs shift parameter-shift -0.642850 0.016461 0.016476 1.0x finite diff h=1e-01 -0.644450 0.141200 0.140847 8.5x finite diff h=1e-02 -0.650000 1.375412 1.371982 83.3x finite diff h=1e-03 -0.315000 13.014324 12.985922 788.2x finite diff h=1e-05 27.500000 1536.488769 1532.901121 93,036.8x

★ Read the $h=10^{-3}$ row. That is the step size the noiseless table called excellent — error $1.07\times10^{-7}$. Under a thousand shots it returns $-0.315$ with a standard deviation of 13.0, against a true value of $-0.644$. The estimator's noise is twenty times the quantity being estimated. It is not a degraded gradient; it is no gradient at all, delivered as a float.

And the ranking inverts. Noiselessly, $h = 10^{-5}$ was the best row in the table and $10^{-1}$ the worst. Under shot noise $10^{-1}$ is the best finite difference and $10^{-5}$ is catastrophic by a factor of 93,000. The sweet spot did not merely worsen — it moved to the opposite end of the table, which is why "tune $h$ on the simulator and deploy" is not a strategy.

The $1/h$ prediction holds:

text predicted std measured std shift 0.017102 0.016461 h=1e-01 0.144345 0.141200 h=1e-03 14.405148 13.014324 h=1e-05 1440.514541 1536.488769

Within about 10% at every row, which for a two-line noise model is as much as one should ask.

Now price the gap. To make a central difference at $h = 10^{-3}$ as precise as the parameter-shift rule already is at 1,000 shots, you would need

text 7.09e+08 shots per evaluation (vs 1,000)

Seven hundred million shots to buy back what $\pi/2$ gives away for free. At $h = 10^{-1}$ it is a milder $7.09\times10^{4}$ — still seventy times the budget, for a method that is also biased by $\sin(h)/h$.

And a warning about the library's own convenience method: diff_method="finite-diff" uses PennyLane's default step and strategy, and under the same 1,000 shots it returned a mean of +18,400 with a standard deviation of 282,000, on a quantity whose true value is $-0.64$. Wrong sign, wrong magnitude, no exception. It is available, it is cheaper in executions, and it is unusable on a sampling device.

The practical rule: if shots is set, the differentiation method is parameter-shift. There is no tuning to do and no trade to consider, and this measurement is why the field settled the question and moved on.

16.4 What a Gradient Costs

The rule needs two evaluations per parameter. That is the entire cost model, and it is worth measuring rather than assuming:

   parameters   circuit executions   per parameter
            1                    3            3.00
            2                    5            2.50
            4                    9            2.25
            8                   17            2.12
           16                   33            2.06

Exactly $2n + 1$ — two shifted evaluations per parameter, plus one forward pass for the value itself.

$$\text{cost of one gradient step} = (2n + 1) \times \text{shots per circuit}$$

This scales linearly in the parameter count, and it is charged on every optimizer iteration. A 100-parameter ansatz optimized for 200 iterations at 4096 shots is $201 \times 200 \times 4096 \approx 1.6 \times 10^8$ shots. Chapter 12 §12.5's point about sessions was about exactly this.

⚙️ Under the Transpiler — PennyLane skips parameters that cannot matter.

Measured on 64 parameters spread across 4 wires, with qml.expval(qml.PauliZ(0)) as the observable:

text 64 params on 4 wires -> 33 executions nonzero gradient entries: 16 of 64

Only the 16 parameters on wire 0 affect $\langle Z_0\rangle$. PennyLane determines this from the circuit structure and evaluates shifts only for those, so the cost is $2 \times 16 + 1 = 33$ rather than $2 \times 64 + 1 = 129$.

A useful optimization, and a useful diagnostic: if your gradient has unexpected zeros, some of your parameters are not connected to your observable. That is usually a bug in the ansatz, not a feature of the landscape — and §16.7 shows a case where it is neither.

💰 Cost and Queue: what $2n+1$ costs in dollars.

§16.4's example is $201 \times 200 \times 4096 = 164{,}659{,}200$ shots. Chapter 39 §39.5 measured two published rates for the same physical work — $0.00035 per shot on superconducting hardware and $0.01 per shot on trapped ions. Apply them:

```text 100 params, 200 iterations, 4,096 shots total shots 164,659,200 of which gradient 163,840,000 = 99.50% of which forward passes 819,200 = 0.50%

at $0.00035/shot   ->  $    57,630.72     (forward passes alone:    $286.72)
at $0.01/shot      ->  $ 1,646,592.00     (forward passes alone:  $8,192.00)

```

★ 99.5% of the bill is the derivative. The thing you actually wanted — the energy — is $287 of a $57,631 invoice. That ratio is $2n/(2n+1)$ and it is above 99% for any $n \geq 100$, so at realistic ansatz sizes "what does this VQE run cost" and "what does its gradient cost" are the same question.

The lever is $n$, and it is a linear lever with an expensive slope:

text 4 params 7,372,800 shots = $ 2,580.48 20 params 33,587,200 shots = $ 11,755.52 100 params 164,659,200 shots = $ 57,630.72

Case Study 2's misdiagnosis is priced here. Going from 4 to 18 parameters to fix what was actually a stopping-criterion bug multiplies the gradient bill by $37/9 \approx 4.1\times$ — every iteration, for the life of the project, in exchange for a slightly worse answer.

And then the queue, which is the larger number. Chapter 39 §39.7 decomposed a 120-iteration VQE run and priced its two execution modes in wall clock: 10.01 hours as 120 separate jobs against 5.52 minutes inside a session, a 108.8× difference for the same circuits and identical shots. That factor is essentially the iteration count, and it stays that way until the device time approaches the queue time. Chapter 12 §12.5 is where the session machinery lives. The shot bill is what you pay; the session is whether you ever see the answer.

One honest caveat, which Chapter 39 §39.5 makes at length: these are list rates with an as_of date, the two providers count differently, and the same computation there priced at $50, $7,432, or $185,542 depending purely on the billing model. Use the arithmetic, re-check the rate.

16.5 An Optimization Loop

The payoff. A two-qubit Hamiltonian:

$$H = Z_0 Z_1 + 0.5\,X_0 + 0.5\,X_1$$

H = qml.Hamiltonian([1.0, 0.5, 0.5],
                    [qml.PauliZ(0) @ qml.PauliZ(1), qml.PauliX(0), qml.PauliX(1)])

@qml.qnode(dev)
def cost(params):
    qml.RY(params[0], 0); qml.RY(params[1], 1)
    qml.CNOT([0, 1])
    qml.RY(params[2], 0); qml.RY(params[3], 1)
    return qml.expval(H)

opt = qml.AdamOptimizer(stepsize=0.1)
for _ in range(400):
    params = opt.step(cost, params)

  exact ground state:  -1.414214
  optimized:           -1.414214      error 6.66e-16

Four parameters, converged to machine precision. The exact answer is $-\sqrt2$, obtained by diagonalizing the Hamiltonian — the reference value Chapter 7 §7.7 insisted on.

🐛 Debug This — "the ansatz is not expressive enough" is usually wrong.

A first attempt at this used GradientDescentOptimizer(stepsize=0.25) for 40 steps and stalled at −1.2485, an error of 0.166 that refused to improve.

The natural diagnosis is that the four-parameter ansatz cannot represent the ground state, and the natural fix is a bigger ansatz. Testing that:

text shallow RY (4 params) -1.414214 error 6.66e-16 RY+RZ (8 params) -1.414212 error 1.33e-06 StronglyEntanglingLayers x3 (18) -1.414213 error 6.32e-08

The four-parameter ansatz was fine. It reaches the exact answer — in fact more precisely than either larger one, which have more directions to wander in. The problem was simply too few steps.

And the instinctive fix is the wrong one. Same ansatz, 200 steps each:

text GradientDescent(0.25) -1.414214 error 2.28e-08 GradientDescent(0.05) -1.259566 error 1.55e-01 Adam(0.1) -1.414214 error 8.76e-09 Momentum(0.1) -1.414214 error 7.80e-10

Reducing the step size makes it worse. The original run at 0.25 was not overshooting — it had not finished. At 0.05 it still has not, after five times as many steps.

Two cheap diagnostics separate a stall from a capacity limit:

Is the trace still descending at the final step? If yes, you ran out of steps. That is a stopping-criterion bug, not a model-capacity problem, and no amount of extra parameters fixes it.

Change the optimizer before changing the ansatz. If a better optimizer closes the gap, the ansatz was never the constraint. Adding parameters to fix an optimizer problem makes it worse — more parameters means a larger gradient bill (§16.4) and, per §16.6, a flatter landscape.

Templates, and where the parameter count comes from

StronglyEntanglingLayers appears in the table above and dominates §16.6, so it is worth knowing what it is. A template is a parameterized sub-circuit with a declared parameter shape, and PennyLane ships several dozen: StronglyEntanglingLayers, BasicEntanglerLayers, AngleEmbedding, AmplitudeEmbedding, QAOAEmbedding, RandomLayers.

The one that matters here is a hardware-efficient ansatz in the standard sense: a layer of arbitrary single-qubit rotations followed by a ring of CNOTs, repeated. Each layer applies Rot(φ, θ, ω) — which PennyLane defines as $R_Z(\omega)\,R_Y(\theta)\,R_Z(\phi)$, three Euler angles spanning all of SU(2) with the rotation gates of Chapter 3 §3.6 — to every wire, then entangles with CNOTs at a per-layer range that cycles.

The parameter shape follows directly, and is worth being able to compute in your head, because §16.4 says the gradient bill is linear in it:

$$\text{params} = 3 \times L \times n$$

Three Euler angles, per wire, per layer. Verified against StronglyEntanglingLayers.shape(6, n):

   wires    shape        params    3 x 6 x n
       2   (6, 2, 3)         36           36
       3   (6, 3, 3)         54           54
       4   (6, 4, 3)         72           72
       6   (6, 6, 3)        108          108
       8   (6, 8, 3)        144          144
      10   (6, 10, 3)       180          180

Those are exactly §16.6's params column. Combine the two facts and the cost model becomes a formula in the ansatz's own hyperparameters: one gradient of a six-layer strongly-entangling ansatz on $n$ wires costs

$$2(3 \times 6 \times n) + 1 = 36n + 1 \text{ circuit executions}$$

At 10 wires that is 361 executions per optimizer step, before any shot count is applied — which is why §16.6's measurement of what those gradients are worth is the load-bearing one.

Templates are also where the surprises live. They are opaque by default — recall §16.1's qml.specs reporting a 54-parameter template as "Total gates: 1" — and they make structural choices on your behalf. The Rot decomposition beginning with $R_Z$ is one such choice, it is entirely reasonable, and §16.7 is what it costs when you measure gradients rather than merely use them. Use templates; know their first and last gate.

16.6 Barren Plateaus

Now the result that governs whether any of this scales.

Set up a standard hardware-efficient ansatz — StronglyEntanglingLayers, six layers — on increasing numbers of qubits. Initialize the parameters randomly, as you would with no better idea, and measure the variance of the gradient across many random initializations.

   qubits   params      Var[grad]    max|grad|
        2       36     1.0431e-01       0.8841
        3       54     5.4316e-02       0.7828
        4       72     2.7964e-02       0.6131
        5       90     1.4716e-02       0.4396
        6      108     8.0048e-03       0.3712
        7      126     3.4924e-03       0.2706
        8      144     1.9115e-03       0.1664
        9      162     9.6043e-04       0.1405
       10      180     4.6517e-04       0.0920

   fit: Var[grad] ~ exp(-0.676 n)   ->  x0.5086 per qubit
   halving every 1.03 qubits
   Var(2 qubits) / Var(10 qubits) = 224x

The gradient variance halves with every qubit added. The fitted factor is 0.5086 per qubit — within 2% of exactly one half, which is the textbook $\mathcal{O}(2^{-n})$ result, measured.

Why this is fatal rather than annoying

A gradient you cannot distinguish from zero is a gradient you cannot follow. And distinguishing a small expectation value from zero costs shots: the statistical error on an expectation value from $N$ shots goes as $1/\sqrt N$, so resolving a gradient of size $g$ needs

$$N \sim 1/g^2$$

Extrapolating the fit:

  at 50 qubits:  Var[grad] ~ 8.76e-16     ->  |grad| ~ 3e-8
  shots required to resolve it:  ~1.1e15

Roughly 10¹⁵ shots — per parameter, per iteration. At any plausible shot rate that is longer than the age of the universe. The optimizer is not slow; it is blind. It sees a flat landscape in every direction and has no information about which way to step.

And note where this bites: 50 qubits is small. It is well below the size at which anyone expects quantum advantage, and comfortably within reach of current hardware. The barren plateau arrives long before the interesting problems do.

Does a shallower ansatz help?

The obvious hope, since deep circuits are more expressive and expressiveness is implicated:

   layers   Var(4 wires)   Var(8 wires)     ratio
        1     1.3842e-01     9.4414e-03     14.7x
        2     4.3994e-02     2.1485e-03     20.5x
        6     2.8436e-02     1.8937e-03     15.0x

No. All three depths fall by a comparable factor over the same four added qubits. Reducing depth raises the overall variance somewhat — a one-layer ansatz starts about 5× higher — but the decay with width is unchanged. You buy a constant factor and the wall stays where it is.

📊 What the Numbers Say: 224× and 88× are the same result, and the difference is the point.

Chapter 32 §32.5 runs this measurement again on a variational classifier and gets a different number — 88× across 2 to 10 qubits, against this chapter's 224×. Two chapters, the same phenomenon, the same qubit range, a 2.5× disagreement. It would be easy to treat one as noise.

It is not noise. It is depth:

text ansatz Var(2q) Var(10q) collapse per qubit Ch.32, StronglyEntangling x2 1.033e-01 1.169e-03 88x 0.571 Ch.16, StronglyEntangling x6 1.0431e-01 4.6517e-04 224x 0.508 Haar-random prediction - - 256x 0.500

Note the first column. At two qubits the two ansätze agree to three digits — 1.033e-01 against 1.0431e-01 — because at two qubits depth barely matters. They diverge as width grows, and they diverge in the direction the theory predicts: the deeper ansatz approaches the asymptotic $\mathcal{O}(2^{-n})$ result (within 13%), and the shallower one falls short of it (by 2.9×) because two layers on ten qubits has not scrambled the state into anything Haar-like yet.

Chapter 32 §32.5 derives this from concentration of measure and is the place to read the argument in full. What it produces is the sentence that governs the field: the ansatz expressive enough to be interesting is the ansatz random enough to be flat. Expressibility and trainability are one knob turned in opposite directions — which is also why the layer sweep above finds a constant-factor improvement and no change in the exponent.

The reading discipline here is the chapter's own. A single barren-plateau number is meaningless without its ansatz and depth. "Gradient variance decays exponentially" is the finding; "×0.51 per qubit" is a property of this circuit, and quoting it as a universal constant would be exactly the error §16.7 catches in a more dramatic form.

🔬 Honest Assessment — what this means for variational quantum algorithms.

VQE and QAOA are the flagship near-term algorithms, and the barren plateau is a measured, reproducible, theoretically-explained obstacle to running them at any interesting scale with random initialization and a hardware-efficient ansatz.

What does not work: picking a shallower ansatz (measured above), more shots (10¹⁵ is not a budget), or a better classical optimizer (the problem is absent information, not poor search).

What is actually proposed, and its status: - Problem-informed ansätze that encode known structure — chemistry-derived ansätze for chemistry problems. Genuinely helps; requires you to already know a lot about your problem. - Smart initialization near a known-good region rather than randomly. Helps, and begs the question of how you found the region. - Layerwise training, growing the circuit gradually. Helps sometimes; has its own plateaus. - Local rather than global observables — measuring $Z_0 Z_1$ rather than a global product. This one has real theoretical backing for shallow circuits, and is why this section's measurement uses a two-qubit observable.

None of these is a general solution, and that is the honest state of the field. Variational algorithms remain worth studying — they are how near-term hardware gets used, and the mitigations above do work on structured problems. But "run VQE on a hardware-efficient ansatz with random initialization at 50 qubits" is not a plan, and the measurement above is why.

Chapter 24 builds VQE properly, with the mitigations, and revisits this honestly.

16.7 Structurally Zero Gradients

A trap discovered while producing §16.6's table, and the reason that table samples all parameters rather than one.

A first attempt measured the gradient of a single fixed parameter — the natural choice — and produced this:

   qubits    Var[grad]     |grad| mean
        2    9.813e-33       7.937e-17
        8    6.942e-34       2.050e-17

Variance of $10^{-33}$ and gradients of $10^{-17}$: machine epsilon. It even decays with qubit count, which looks like an extremely convincing barren plateau.

It is not a barren plateau. It is one gate.

StronglyEntanglingLayers begins with Rot(φ, θ, ω) on each wire, and Rot starts with an RZ. The first layer acts on $|0\rangle$, and

$$R_Z(\phi)|0\rangle = e^{-i\phi/2}|0\rangle$$

a global phase, which is unobservable. The derivative with respect to that parameter is not small; it is identically zero, for every input, at every qubit count.

Confirmed directly:

  g[0,0,0] (first RZ of first Rot, on |0>)  = 9.714e-17     <- structurally zero
  g[0,:,0] (that angle on every wire)       = [0. 0. 0.]
  g[2,1,1] (a middle-layer angle)           = -0.143103     <- a real gradient
  max |gradient| over all 36 params         =  0.495574

The real gradients are $\mathcal{O}(0.1)$ — fourteen orders of magnitude larger than the artifact.

🐛 Debug This — a zero gradient has three possible causes, and they need different responses.

Cause	Signature	Response
Structurally zero	exactly ~1e-17, at every parameter value and every size	fix the ansatz, or exclude the parameter
Disconnected from the observable	exactly zero, and the parameter's gates are on unmeasured wires	ansatz bug (§16.4)
Barren plateau	small but nonzero, shrinking exponentially with width	the real problem

The discriminator is cheap: change one qubit and one parameter value. A structurally zero gradient stays ~1e-17 regardless. A barren plateau's gradient is a genuine random variable that happens to be small — nonzero, differently sized on each draw, and shrinking with $n$.

Had this gone uncaught, the chapter would have reported a barren plateau eighteen orders of magnitude more severe than the real one, with a clean exponential fit to make it convincing.

Sample all parameters, not one. And treat any measured value near machine epsilon as an artifact until proven otherwise — this book's recurring lesson that a number can be precise, reproducible, and about something other than what you think.

How many parameters are dead? A census

The first RZ is the famous one, and it is not the only one. Since the correct measurement filters structurally-zero entries anyway, it costs nothing to ask how many there are — and the answer is larger, more structured, and less predictable than the single-gate story suggests.

Take the same setup as §16.6 — six layers, observable $Z_0 \otimes Z_1$ — and mark every parameter whose gradient is $\leq 10^{-14}$ on every draw:

   wires   params   first layer   last layer   elsewhere   total    % dead
       2       36             2            4           0       6     16.7%
       3       54             3            7           0      10     18.5%
       4       72             4            8           6      18     25.0%
       5       90             5            9           0      14     15.6%
       6      108             6            8           0      14     13.0%
       8      144             8           16           1      25     17.4%
      10      180            10           22          13      45     25.0%

Between 13% and 25% of the parameters are structurally dead, and the counts are identical across 8 draws and 20 draws, at seeds 11 and 999. These are not small gradients. They are zeros.

The first-layer column is the known effect and behaves exactly as §16.7 says: exactly $n$ of them, one per wire, always parameter index 0 — the leading $R_Z$ acting on $|0\rangle$.

The last-layer column is a second, mirror-image effect, and it is derivable. The final Rot layer is followed only by CNOTs before a $Z$-basis expectation value. Conjugating by a CNOT maps $Z$-type Paulis to $Z$-type Paulis — $Z_c \to Z_c$ and $Z_t \to Z_c Z_t$ — so the observable propagated back to just after the last rotations is still a product of $Z$s. And $R_Z$ commutes with any product of $Z$s. Therefore the trailing $R_Z$ of the last Rot (parameter index 2) has identically zero derivative on every wire, for the same reason the leading one does: at the input the state is a $Z$ eigenstate; at the output the observable is a $Z$ operator.

Propagating $Z_0 \otimes Z_1$ back through the last layer's CNOT ring gives more than that:

   wires   last-layer range   effective observable before the ring
       2         r=1          Z0
       3         r=2          Z2
       4         r=3          Z0 (x) Z2
       6         r=1          Z0 (x) Z2 (x) Z3 (x) Z4 (x) Z5

Wires missing from that list carry an identity, and a rotation against an identity observable has zero gradient in all three angles, not just the trailing one. That is a §16.4-style disconnection — produced not by a buggy ansatz but by the standard template's own entangling ring. Predicting the last-layer mask this way and comparing to the measurement:

   wires    predicted zeros   measured zeros   match
       2                  4                4   yes
       3                  7                7   yes
       4                  8                8   yes
       6                  8                8   yes

Exact, at every size. Two of §16.7's three categories turn out to be the same mechanism seen from opposite ends of the circuit.

Now the useful part, which is the elsewhere column. It is 0, 0, 6, 0, 0, 1, 13 — not monotonic, not proportional to width, and not predictable from the ansatz's name. At four wires the range-3 ring happens to split into pairs that leave two wires disconnected a full layer earlier; at ten wires thirteen parameters die somewhere in the middle. There is no formula. There is a measurement.

⚠️ Common Pitfall: reporting a gradient variance without reporting the denominator.

A quarter of the parameters being dead does not change §16.6's variance — the filter removes them — but it changes anything computed per parameter. §16.4's gradient bill is $2n+1$ where $n$ is the number of trainable parameters, and at 10 wires 45 of the 180 contribute nothing: you are paying for 361 circuit executions and 90 of them measure a quantity that is zero by construction.

That is 25% of the shot budget spent on arithmetic whose answer is known in advance, which at §16.4's rates is $\$14{,}400$ of a $\$57{,}631$ run.

The two failure modes point in opposite directions and both are easy. Silently dropping the zeros inflates any "average gradient magnitude" you report. Silently keeping them deflates it, and deflates it more at larger widths, since the dead fraction is not constant — which would manufacture a barren plateau out of bookkeeping. Case Study 1 is the extreme version of the second error.

measure_gradient_variance in the project code reports the exclusion count for exactly this reason. A filter that announces what it dropped invites the question of what those values were; a silent one does not.

🧱 Project Checkpoint — vqelab/variational.py: the optimization layer.

vqelab can choose good qubits (Ch. 12), mitigate what is left (Ch. 13), cross framework boundaries (Ch. 14), and price the result fault-tolerantly (Ch. 15). This adds the loop that actually finds an answer.

gradient_cost(n_params, shots, iterations) returns the total shot count from §16.4's $2n+1$ rule, so the bill is visible before you commit to it rather than after.

measure_gradient_variance(ansatz, n_wires, samples) reproduces §16.6's measurement for your ansatz, excluding structurally-zero parameters — with the exclusion logged, not silent, because §16.7 is exactly how a silent exclusion would have lied.

classify_zero_gradient(qnode, params, index) implements §16.7's three-way discriminator, returning "structural", "disconnected", or "plateau" with the evidence for the call.

optimize(cost, params, optimizer, ...) runs the loop and returns a VariationalResult carrying the energy trace, the shot cost, and — critically — a converged flag that is false when the trace is still descending at the final step, because §16.5's stall was a stopping-criterion bug misdiagnosed as an ansatz problem.

Its tests assert the two measured facts: a gradient costs exactly $2n+1$ executions, and gradient variance decays by more than 100× from 2 to 10 qubits.

16.8 Where PennyLane Fits

PennyLane's plugin system is worth noting: pennylane-qiskit and pennylane-cirq let a QNode execute on those backends, so PennyLane is often used as a differentiation layer on top of another framework rather than as a replacement for one. That is arguably its best use — Chapter 18 returns to it.

Broadcasting: the throughput half of the framework

The parameter-shift rule is why PennyLane exists. Broadcasting is why anything written in it finishes.

A QNode accepts a leading batch axis on its arguments and evaluates the whole array in one pass. Measured on 256 samples through AngleEmbedding + CNOT, returning $\langle Z_0\rangle$:

   256 samples, looped     256 executions      608.15 ms
   256 samples, batched    256 executions        2.86 ms      212x wall clock
   identical results:  max absolute difference  0.00e+00

★ Read the middle column, not the right one. The execution count is the same: 256 either way. Broadcasting did not make the quantum work cheaper — it removed about 2.4 milliseconds of Python and PennyLane dispatch per call, paid 256 times. On a two-qubit simulator the dispatch dominates the circuit by orders of magnitude, so removing it looks like a 212× speedup, and it is one, and it is not a quantum effect at all. (The exact factor is machine-load dependent; Chapter 34 §34.2 measured 478× and 868× on repeat runs of the same experiment.)

This is one of the book's most-repeated measurements, in four settings. Chapter 33 §33.3 got ~100× by batching a training loop; Chapter 34 §34.2 got 478× on a 40,401-entry Gram matrix, from 1,270 µs per entry to 2.7 µs; Chapter 39 §39.3 got ~99× on hardware by submitting 100 circuits as one job instead of 100. The loop you wrote is not the cost you think it is.

And here is the part that matters for §16.4's bill. Take the gradient of a batched loss — 3 parameters, 256 samples:

   1,792 executions  =  7 x 256  =  (2n+1) x batch size

Broadcasting does not amortize the parameter-shift cost across the batch. Each sample needs its own pair of shifted evaluations, because each sample is a different circuit. Batching removes the per-submission cost and leaves the per-shot cost exactly where it was — which is Chapter 34 §34.2's conclusion arriving here by a different route. Batching makes the gradient possible; it does not make it cheap.

Interfaces, and when the choice matters

qml.qnode takes an interface argument: "autograd" (the default, and what pennylane.numpy provides), "torch", "jax", "tensorflow", "numpy", or "auto". In this book's environment — Appendix C ships none of PyTorch, JAX, or TensorFlow — only autograd, numpy, and auto resolve; the others raise on import. No timing comparison between interfaces is offered here, because none was run.

What the setting actually does is worth stating precisely, because it is routinely misunderstood. The interface decides which framework owns the autodiff tape, not how the derivative is computed. Choosing "torch" makes the QNode a node in a PyTorch graph, so gradients flow through it into classical layers on either side and loss.backward() reaches your circuit parameters. It does not change diff_method. On a sampling device the derivative is still parameter-shift and still costs $2n+1$ circuit executions per gradient. An interface cannot make a quantum gradient cheaper; only a smaller $n$ can.

So the choice is a systems question rather than a performance one:

• The circuit is the whole model. Autograd and qml.AdamOptimizer are sufficient, and this chapter's §16.5 is the shape of it. Adding torch buys a dependency and a device-placement question for nothing.
• The circuit is one layer inside a classical network. The interface is the entire point — it is what makes end-to-end training work rather than requiring you to hand-wire a chain rule across the boundary.
• The classical side is the expensive side. JAX's jit and vmap are worth reaching for, on the classical half.

Which is a decision about where the hybrid boundary sits, and Chapter 35 §35.6 is the chapter that takes that question seriously — it identifies three places the boundary can go and argues that the choice is usually made by accident.

🔀 In Another Framework: nobody else ships the gradient.

Chapter 24 §24.2 compares the full optimization loop across frameworks and finds Qiskit + COBYLA and PennyLane + Adam both landing on $-1.857275030$ Ha. The interesting comparison here is narrower: who gives you $\partial f/\partial\theta$?

Qiskit 2.5.1 — nothing in core. Searching the top-level namespace for anything gradient-related returns an empty list. You get Parameter, assign_parameters, and the Estimator primitives, and then you write the parameter-shift bookkeeping yourself or install the separate qiskit-algorithms package for its gradient classes. Neither qiskit-algorithms nor qiskit-machine-learning is in the metapackage — pip install qiskit does not get you a derivative. This is why Chapter 24's Qiskit path uses COBYLA, which is gradient-free: not because gradient-free is better, but because it is what the default install supports.

Cirq 1.7.0 — nothing, and deliberately so. PhaseGradientGate is a gate, not a differentiation API. Cirq has PauliSum and a fast simulator and stops there, which is Chapter 14 §14.1's trade stated exactly: the framework that gives you the least gives you the fewest surprises about what is happening.

Q# (Chapter 15) — no autodiff and no classical optimizer in the language. Q# is aimed at resource estimation and fault-tolerant reasoning, and a variational loop is not what it is for.

This is the strongest single argument for PennyLane in the book. The other three frameworks each do something PennyLane does not, and Chapter 18 is about combining them. But if you need a gradient on hardware, exactly one of the four hands it to you — and via pennylane-qiskit it will hand it to you while running on Qiskit's backends, which is the arrangement §16.8's table is really recommending.

🧪 Run It: three experiments, twenty minutes, no hardware.

1. Break the gradient on purpose. Take §16.1's QNode and call qml.grad on np.array(0.7) instead of pnp.array(0.7, requires_grad=True). Confirm you get () and not an exception. Then set warnings.simplefilter("error", UserWarning) and watch the same call become a traceback. That two-line habit is worth more than anything else in this chapter on your first real project.

2. Catch the framework choosing for you. Build the same QNode twice, once with shots=None and once with @qml.set_shots(1000), and print qml.workflow.get_best_diff_method(qnode)(theta) for each. You should see backprop and parameter-shift. Then ask for diff_method="backprop" with shots set and read the error message carefully — it is the clearest one-sentence explanation of why hardware gradients are expensive that any framework produces.

3. Watch the sweet spot move. Reproduce §16.3's finite-difference table on an exact device, note which $h$ wins. Then re-run it under @qml.set_shots(1000), repeated 200 times per step size, and report the standard deviation rather than a single value. The winner changes ends of the table. This is the noise report above, and running it yourself is the fastest way to internalise why the parameter-shift rule is not merely a tidier finite difference.

All three run on default.qubit in seconds. The third one is also Exercise 16.8.

16.9 Summary

A QNode is a circuit bound to a device, and it is differentiable. The function body is the circuit; the return declares the measurement; pennylane.numpy tracks gradients. There is no circuit object and no Parameter class — so Chapter 8's parameter-ordering bug cannot occur here either, for a third distinct reason.

★ The parameter-shift rule gives exact gradients from circuit evaluations:

$$\frac{\partial f}{\partial\theta} = \frac{f(\theta + \pi/2) - f(\theta - \pi/2)}{2}$$

Verified to zero error against $-\sin(0.7)$. Exact because the expectation value is exactly sinusoidal in a Pauli-generated rotation angle — not an approximation, and with no step size to tune. Finite differences by contrast have a sweet spot and degrade on both sides of it (error $1.07\times10^{-3}$ at $h{=}10^{-1}$, $1.98\times10^{-11}$ at $10^{-5}$, then worse again at $10^{-7}$).

A gradient costs exactly $2n+1$ circuit executions, measured: 16 parameters → 33 executions. Charged on every optimizer iteration — a 100-parameter, 200-iteration, 4096-shot run is $1.6\times10^8$ shots. PennyLane skips parameters that cannot affect the observable (64 params across 4 wires → 33 executions, 16 nonzero gradients), which is also a useful ansatz bug detector.

Diagnose optimizer failures before enlarging the ansatz. A four-parameter ansatz stalled at −1.2485 under GradientDescentOptimizer for 40 steps; the same ansatz under Adam for 400 steps reached −1.414214, error 6.7e−16 — machine precision, and better than the 8- and 18-parameter alternatives. Change the optimizer and the step count first; adding parameters to fix an optimizer problem makes both the gradient bill and the landscape worse.

★★ Barren plateaus are real, and measured. Gradient variance for a randomly-initialized hardware-efficient ansatz decays as $\exp(-0.676\,n)$ — ×0.51 per qubit, halving every 1.03 qubits, a 224× drop from 2 to 10 qubits. That is the textbook $\mathcal{O}(2^{-n})$ result.

Since resolving a gradient of size $g$ needs $N \sim 1/g^2$ shots, at 50 qubits the requirement is about 10¹⁵ shots per parameter per iteration. The optimizer is not slow — it is blind. And 50 qubits is small: the plateau arrives well before the problems worth solving.

Shallower ansätze do not escape it. One, two, and six layers all fall by a comparable factor over the same four added qubits (14.7×, 20.5×, 15.0× from 4 to 8 qubits). You buy a constant factor in the overall scale; the decay with width does not move.

★ Distinguish a barren plateau from a structurally zero gradient. StronglyEntanglingLayers opens with RZ acting on $|0\rangle$ — a global phase — so that parameter's derivative is identically zero, ~1e-17, at every size. Sampling it produced a beautifully exponential fit to an artifact eighteen orders of magnitude more severe than the real effect. Sample all parameters, and treat any value near machine epsilon as an artifact until proven otherwise.

PennyLane is the framework for variational and machine-learning work, and its plugin system makes it most valuable as a differentiation layer on top of Qiskit or Cirq rather than a replacement.

Next: Chapter 17 — Amazon Braket, which takes a different position again: not a better abstraction over quantum circuits, but a single API over other people's hardware, including trapped ions and neutral atoms whose gate sets and connectivity differ from everything in Part II.