Gradients

Six estimators behind one grad(), plus the shift-rule machinery they share.

qmlkit.gradients.batch

Gradients for a whole training batch in one pass. param_shift_grad_batch never inspects a state, so it works on every backend including a sampling-only device, and is the batched submission real hardware wants.

batch

Gradients for a whole training batch, in one pass rather than one per sample.

A training step evaluates the same circuit structure at one parameter vector per sample. The forward pass has been batched for a while; the backward pass had not, and the backward pass is where the time goes — measured on a 4-qubit VQC at batch 128, the forward pass was 13 ms of an 1170 ms step. 99% of training was still running one sample at a time.

Two routes, and the difference between them is the point of this module.

:func:param_shift_grad_batch is backend-agnostic. A shift rule only ever needs the circuit run at shifted angles, so a whole batch's gradient is one big set of evaluations — batch x 2P of them — with no inspection of the state at all. It goes through :meth:~qmlkit.core.backends.base.Backend.expectation_over_slots, which every backend has, so this works on NumPy, Qiskit, Cirq, SpinQit, a sampling-only device, and anything registered later. On hardware it is exactly the batched submission a provider wants: one job instead of batch x 2P blocking calls.

:func:adjoint_grad_batch is faster and simulator-only, because it reads the state.

Both are exact, and both are asserted equal to the per-sample functions they replace.

The shifts happen in slot space, not logical-parameter space, which is why the batch primitives underneath are slot-based. A shift rule moves one occurrence of a parameter and a weight-tied parameter has several; a batched routine written against logical parameters cannot express that, and would silently compute a different derivative for any circuit with shared weights.

param_shift_grad_batch

param_shift_grad_batch(
    spec: CircuitSpec,
    thetas: ArrayLike,
    obs: Observable | None = None,
    shots: int | None = None,
    backend: Backend | str | None = None,
    seed: int | None = None,
) -> NDArray[Any]

Parameter-shift gradients for a batch, as one set of evaluations.

Returns (batch, n_params).

Every shifted circuit the whole batch needs is assembled first and handed to the backend in a single call, so a backend that can evaluate many circuits at once — or submit them as one job — gets to. Nothing here inspects a state, so this is valid on hardware and under sampling.

Weight-tied parameters accumulate across occurrences and rescaled references get their chain-rule factor, exactly as the per-sample :func:~qmlkit.gradients.parameter_shift.param_shift_grad does.

Source code in src/qmlkit/gradients/batch.py

def param_shift_grad_batch(
    spec: CircuitSpec,
    thetas: ArrayLike,
    obs: Observable | None = None,
    shots: int | None = None,
    backend: Backend | str | None = None,
    seed: int | None = None,
) -> npt.NDArray[Any]:
    """Parameter-shift gradients for a batch, as **one** set of evaluations.

    Returns ``(batch, n_params)``.

    Every shifted circuit the whole batch needs is assembled first and handed to the
    backend in a single call, so a backend that can evaluate many circuits at once —
    or submit them as one job — gets to. Nothing here inspects a state, so this is
    valid on hardware and under sampling.

    Weight-tied parameters accumulate across occurrences and rescaled references get
    their chain-rule factor, exactly as the per-sample
    :func:`~qmlkit.gradients.parameter_shift.param_shift_grad` does.
    """
    obs = Z(0) if obs is None else obs
    be = get_backend(backend)
    angles = spec.bind_slots_batch(np.atleast_2d(np.asarray(thetas, dtype=float)))
    batch = angles.shape[0]
    slots = spec.slots()
    if not slots:
        return np.zeros((batch, spec.n_params), dtype=float)

    # ---- assemble every shifted row the batch needs, once ------------------ #
    blocks: list[npt.NDArray[Any]] = []
    plan: list[tuple[int, float]] = []  # (slot index, coefficient)
    needs_unshifted = False
    for i, slot in enumerate(slots):
        rule = rule_for_gate(slot.gate)
        for shift, coeff in zip(rule.shifts, rule.coeffs, strict=False):
            shifted = angles.copy()
            shifted[:, i] += shift
            blocks.append(shifted)
            plan.append((i, coeff))
        needs_unshifted = needs_unshifted or rule.needs_unshifted
    if needs_unshifted:
        blocks.append(angles)

    values = be.expectation_over_slots(
        spec, np.concatenate(blocks, axis=0), obs, shots, seed
    ).reshape(len(blocks), batch)

    # ---- fold the evaluations back into per-parameter gradients ------------ #
    grad = np.zeros((batch, spec.n_params), dtype=float)
    for k, (i, coeff) in enumerate(plan):
        ref = slots[i].ref
        grad[:, ref.index] += coeff * values[k] * ref.scale
    if needs_unshifted:
        base = values[-1]
        for slot in slots:
            rule = rule_for_gate(slot.gate)
            if rule.needs_unshifted:
                grad[:, slot.ref.index] += rule.unshifted_coeff * base * slot.ref.scale
    return grad

param_shift_batch_cost

param_shift_batch_cost(
    spec: CircuitSpec, batch: int
) -> int

Circuit evaluations a batched parameter-shift gradient submits, in total.

Source code in src/qmlkit/gradients/batch.py

def param_shift_batch_cost(spec: CircuitSpec, batch: int) -> int:
    """Circuit evaluations a batched parameter-shift gradient submits, in total."""
    per_sample = sum(rule_for_gate(s.gate).n_evaluations for s in spec.slots())
    return per_sample * batch

adjoint_grad_batch

adjoint_grad_batch(
    spec: CircuitSpec,
    thetas: ArrayLike,
    obs: Observable | None = None,
    backend: Backend | str | None = None,
) -> NDArray[Any]

Adjoint gradients for a batch in one forward and one backward sweep.

Returns (batch, n_params). Exact, and independent of P in cost, but it needs the statevector — so it is simulator-only, and :func:param_shift_grad_batch is what a device uses.

The sweep is the one in :mod:~qmlkit.gradients.adjoint with the batch carried as a leading axis: every gate is undone for the whole stack at once, and each parameter's contribution comes out as one inner product per row.

Source code in src/qmlkit/gradients/batch.py

def adjoint_grad_batch(
    spec: CircuitSpec,
    thetas: ArrayLike,
    obs: Observable | None = None,
    backend: Backend | str | None = None,
) -> npt.NDArray[Any]:
    """Adjoint gradients for a batch in one forward and one backward sweep.

    Returns ``(batch, n_params)``. Exact, and independent of ``P`` in cost, but it
    needs the statevector — so it is simulator-only, and
    :func:`param_shift_grad_batch` is what a device uses.

    The sweep is the one in :mod:`~qmlkit.gradients.adjoint` with the batch carried as
    a leading axis: every gate is undone for the whole stack at once, and each
    parameter's contribution comes out as one inner product per row.
    """
    obs = Z(0) if obs is None else obs
    be = get_backend(backend)
    if not be.supports_statevector:
        raise ValueError(
            f"the {be.name!r} backend has no statevector, so it cannot differentiate by "
            'the adjoint method; use method="parameter-shift", which batches too'
        )

    values = np.atleast_2d(np.asarray(thetas, dtype=float))
    angles = spec.bind_slots_batch(values)
    batch, n = angles.shape[0], spec.n_qubits
    slots = spec.slots()

    psi = be.statevector_batch_slots(spec, angles).reshape((batch,) + (2,) * n)
    lam = _apply_observable_batch(psi, obs, n)
    grad = np.zeros((batch, spec.n_params), dtype=float)

    # walk the slots backwards in step with the ops, the way the scalar version does
    cursor = len(slots)
    for op_index in range(len(spec.ops) - 1, -1, -1):
        op = spec.ops[op_index]
        first = cursor
        while first > 0 and slots[first - 1].op_index == op_index:
            first -= 1
        columns = angles[:, first:cursor] if first < cursor else None
        cursor = first

        u = _batched_matrices(op, columns)
        psi = _apply_batch(psi, _dagger(u), op.qubits)  # the state before this gate

        if columns is not None and columns.shape[1] == 1:
            du = _batched_derivatives(op, columns)
            mu = _apply_batch(psi, du, op.qubits)
            overlap = np.einsum("bi,bi->b", lam.reshape(batch, -1).conj(), mu.reshape(batch, -1))
            ref = slots[first].ref
            grad[:, ref.index] += 2.0 * np.real(overlap) * ref.scale
        elif columns is not None and columns.shape[1] > 1:  # pragma: no cover
            raise NotImplementedError(
                f"gate {op.gate!r} has {columns.shape[1]} parameters; batched adjoint "
                'handles one per gate. Use method="parameter-shift".'
            )

        lam = _apply_batch(lam, _dagger(u), op.qubits)
    return grad

grad_batch

grad_batch(
    spec: CircuitSpec,
    thetas: ArrayLike,
    obs: Observable | None = None,
    method: str = "auto",
    shots: int | None = None,
    backend: Backend | str | None = None,
    seed: int | None = None,
) -> NDArray[Any]

(batch, n_params) gradients, by whichever route fits the backend.

"auto" picks adjoint on an exact simulator and parameter-shift otherwise — the same rule :func:~qmlkit.gradients.dispatch.choose_method uses for one sample, since sampling puts the statevector out of reach by definition.

>>> import numpy as np, qmlkit as qk
>>> a = qk.hardware_efficient(3, 2)
>>> qk.grad_batch(a.build(), np.zeros((5, a.n_params)), qk.Z(0)).shape
(5, 12)

Source code in src/qmlkit/gradients/batch.py

def grad_batch(
    spec: CircuitSpec,
    thetas: ArrayLike,
    obs: Observable | None = None,
    method: str = "auto",
    shots: int | None = None,
    backend: Backend | str | None = None,
    seed: int | None = None,
) -> npt.NDArray[Any]:
    """``(batch, n_params)`` gradients, by whichever route fits the backend.

    ``"auto"`` picks adjoint on an exact simulator and parameter-shift otherwise —
    the same rule :func:`~qmlkit.gradients.dispatch.choose_method` uses for one
    sample, since sampling puts the statevector out of reach by definition.

        >>> import numpy as np, qmlkit as qk
        >>> a = qk.hardware_efficient(3, 2)
        >>> qk.grad_batch(a.build(), np.zeros((5, a.n_params)), qk.Z(0)).shape
        (5, 12)
    """
    if method == "auto":
        from qmlkit.gradients.adjoint import supports_adjoint

        method = (
            "adjoint" if shots is None and supports_adjoint(spec, backend) else "parameter-shift"
        )
    if method == "adjoint":
        if shots is not None:
            raise ValueError("adjoint differentiation is exact; it cannot take shots")
        return adjoint_grad_batch(spec, thetas, obs, backend)
    if method == "parameter-shift":
        return param_shift_grad_batch(spec, thetas, obs, shots, backend, seed)

    from qmlkit.utils.errors import unknown

    raise unknown(
        "batched gradient method",
        method,
        ("auto", "adjoint", "parameter-shift"),
        hint=" Other methods have no batched form yet; loop over qk.grad for those.",
    )

qmlkit.gradients.dispatch

dispatch

One gradient function, several methods, and a registry for adding more.

qk.grad(spec, theta, obs)                      # picks the right method for you
qk.grad(spec, theta, obs, method="parameter-shift")

method="auto" (the default) uses adjoint when every gate has a closed-form derivative and the backend can hand back a statevector — exact, and independent of the parameter count. It falls back to parameter-shift otherwise, which is exact too, just 2P circuits instead of one pass.

A researcher with a different estimator registers it and it becomes a keyword everywhere the library takes method=:

@register_gradient("my_estimator")
def my_estimator(spec, theta, obs, *, backend=None, shots=None, **kw):
    ...

register_gradient

register_gradient(
    name: str, fn: GradFn | None = None
) -> Callable[[GradFn], GradFn] | GradFn

Register a gradient estimator under name. Usable as a decorator.

Source code in src/qmlkit/gradients/dispatch.py

def register_gradient(name: str, fn: GradFn | None = None) -> Callable[[GradFn], GradFn] | GradFn:
    """Register a gradient estimator under ``name``. Usable as a decorator."""

    def _register(f: GradFn) -> GradFn:
        if name in _METHODS:
            raise ValueError(f"gradient method {name!r} is already registered")
        _METHODS[name] = f
        return f

    return _register if fn is None else _register(fn)

choose_method

choose_method(
    spec: CircuitSpec,
    backend: Backend | str | None = None,
    shots: int | None = None,
) -> str

What method="auto" resolves to, and why.

Sampling means the statevector is off the table by definition, so shots force parameter-shift.

Source code in src/qmlkit/gradients/dispatch.py

def choose_method(
    spec: CircuitSpec, backend: Backend | str | None = None, shots: int | None = None
) -> str:
    """What ``method="auto"`` resolves to, and why.

    Sampling means the statevector is off the table by definition, so shots force
    parameter-shift.
    """
    if shots is not None:
        return "parameter-shift"
    from qmlkit.gradients.adjoint import supports_adjoint

    return "adjoint" if supports_adjoint(spec, backend) else "parameter-shift"

grad

grad(
    spec: CircuitSpec,
    theta: ArrayLike,
    obs: Observable | None = None,
    method: str = "auto",
    backend: Backend | str | None = None,
    shots: int | None = None,
    **kwargs: object,
) -> NDArray[Any]

Gradient of <obs> with respect to theta.

Parameters:

Name	Type	Description	Default
method	str	auto (default), adjoint, parameter-shift, spsa, or finite-diff. Registered methods are accepted by name too.	'auto'
shots	int | None	Sampling budget. Anything other than None rules out adjoint.	None

Source code in src/qmlkit/gradients/dispatch.py

def grad(
    spec: CircuitSpec,
    theta: ArrayLike,
    obs: Observable | None = None,
    method: str = "auto",
    backend: Backend | str | None = None,
    shots: int | None = None,
    **kwargs: object,
) -> npt.NDArray[Any]:
    """Gradient of ``<obs>`` with respect to ``theta``.

    Parameters
    ----------
    method
        ``auto`` (default), ``adjoint``, ``parameter-shift``, ``spsa``, or
        ``finite-diff``. Registered methods are accepted by name too.
    shots
        Sampling budget. Anything other than ``None`` rules out adjoint.
    """
    obs = Z(0) if obs is None else obs
    resolved = choose_method(spec, backend, shots) if method == "auto" else method
    try:
        fn = _METHODS[resolved]
    except KeyError:
        raise unknown(
            "gradient method",
            resolved,
            list_gradient_methods(),
            hint='"auto" chooses for you; add your own with register_gradient(name, fn).',
            error=KeyError,
        ) from None
    return fn(spec, np.asarray(theta, dtype=float), obs, backend=backend, shots=shots, **kwargs)

hessian

hessian(
    spec: CircuitSpec,
    theta: Sequence[float],
    obs: Observable | None = None,
    backend: Backend | str | None = None,
    eps: float = 0.0001,
) -> NDArray[Any]

Second derivatives, by differencing the exact gradient.

The gradient itself is exact (adjoint), so only the outer derivative is approximated — far more accurate than differencing the expectation twice.

Source code in src/qmlkit/gradients/dispatch.py

def hessian(
    spec: CircuitSpec,
    theta: Sequence[float],
    obs: Observable | None = None,
    backend: Backend | str | None = None,
    eps: float = 1e-4,
) -> npt.NDArray[Any]:
    """Second derivatives, by differencing the exact gradient.

    The gradient itself is exact (adjoint), so only the outer derivative is
    approximated — far more accurate than differencing the expectation twice.
    """
    obs = Z(0) if obs is None else obs
    arr = np.asarray(theta, dtype=float).ravel()
    p = arr.size
    out = np.zeros((p, p))
    for k in range(p):
        plus, minus = arr.copy(), arr.copy()
        plus[k] += eps
        minus[k] -= eps
        out[k] = (
            grad(spec, plus, obs, backend=backend) - grad(spec, minus, obs, backend=backend)
        ) / (2 * eps)
    return 0.5 * (out + out.T)  # symmetrise away the differencing asymmetry

gradient_cost

gradient_cost(
    spec: CircuitSpec, method: str = "parameter-shift"
) -> int | str

Circuit evaluations one gradient needs under a given method.

Source code in src/qmlkit/gradients/dispatch.py

def gradient_cost(spec: CircuitSpec, method: str = "parameter-shift") -> int | str:
    """Circuit evaluations one gradient needs under a given method."""
    from qmlkit.gradients.hadamard import hadamard_grad_cost
    from qmlkit.gradients.parameter_shift import grad_circuit_cost

    costs = {
        "adjoint": 1,
        "backprop": 1,
        "hadamard": hadamard_grad_cost(spec),
        "parameter-shift": grad_circuit_cost(spec),
        "finite-diff": 2 * spec.n_params,
        "spsa": 2,
    }
    return costs.get(method, "unknown")

qmlkit.gradients.rules

rules

Shift rules, derived rather than remembered.

A circuit expectation as a function of one gate angle is a finite Fourier series whose frequencies are the unique positive gaps between the eigenvalues of that gate's generator:

.. math:: f(\theta) = a_0 + \sum_{k} a_k \cos(\omega_k \theta) + b_k \sin(\omega_k \theta)

We want coefficients :math:c_i and shifts :math:s_i with :math:\sum_i c_i f(\theta + s_i) = f'(\theta) for every such :math:f. Expanding and matching the :math:\cos(\omega\theta) and :math:\sin(\omega\theta) terms gives, for each frequency :math:\omega:

.. math:: \sum_i c_i \cos(\omega s_i) = 0, \qquad \sum_i c_i \sin(\omega s_i) = \omega

Choosing antisymmetric shifts :math:\{+s_1, -s_1, \dots\} with antisymmetric coefficients satisfies the cosine equation and :math:\sum_i c_i = 0 identically, leaving an :math:R \times R linear system in the positive half. We solve it.

Deriving the rule instead of hardcoding constants means a new gate needs only its frequencies declared — no new gradient code, and no chance of a transcribed constant being subtly wrong. :func:four_term_rule reproduces the textbook controlled-rotation constants, and the test suite asserts exactly that.

ShiftRule dataclass

ShiftRule(
    shifts: tuple[float, ...],
    coeffs: tuple[float, ...],
    unshifted_coeff: float = 0.0,
)

f'(theta) = sum_i coeffs[i] * f(theta + shifts[i]) (+ an unshifted term).

general_shift_rule

general_shift_rule(
    frequencies: Sequence[float],
    shifts: Sequence[float] | None = None,
) -> ShiftRule

Build the exact shift rule for a generator with these frequencies.

Source code in src/qmlkit/gradients/rules.py

def general_shift_rule(
    frequencies: Sequence[float], shifts: Sequence[float] | None = None
) -> ShiftRule:
    """Build the exact shift rule for a generator with these frequencies."""
    freqs = tuple(sorted({float(f) for f in frequencies}))
    if not freqs:
        raise ValueError("a shift rule needs at least one generator frequency")
    if any(f <= 0 for f in freqs):
        raise ValueError(f"frequencies must be positive, got {freqs}")

    pos = np.asarray(shifts, dtype=float) if shifts is not None else _positive_shifts(freqs)
    if pos.size != len(freqs):
        raise ValueError(f"need exactly {len(freqs)} positive shifts, got {pos.size}")

    # 2 * sum_j c_j sin(w_k s_j) = w_k
    a = 2.0 * np.sin(np.outer(np.asarray(freqs), pos))
    if abs(np.linalg.det(a)) < 1e-12:
        raise ValueError(
            f"degenerate shift choice {pos.tolist()} for frequencies {list(freqs)}; "
            "pick different shifts"
        )
    c = np.linalg.solve(a, np.asarray(freqs, dtype=float))

    out_shifts: list[float] = []
    out_coeffs: list[float] = []
    for s, ci in zip(pos, c, strict=False):
        out_shifts += [float(s), float(-s)]
        out_coeffs += [float(ci), float(-ci)]
    return ShiftRule(tuple(out_shifts), tuple(out_coeffs))

rule_for_frequencies cached

rule_for_frequencies(
    frequencies: tuple[float, ...],
) -> ShiftRule

Cached :func:general_shift_rule — rules are pure functions of the spectrum.

Source code in src/qmlkit/gradients/rules.py

@cache
def rule_for_frequencies(frequencies: tuple[float, ...]) -> ShiftRule:
    """Cached :func:`general_shift_rule` — rules are pure functions of the spectrum."""
    return general_shift_rule(frequencies)

two_term_rule

two_term_rule() -> ShiftRule

The familiar Pauli-rotation rule: shifts +-pi/2, coefficients +-1/2.

Source code in src/qmlkit/gradients/rules.py

def two_term_rule() -> ShiftRule:
    """The familiar Pauli-rotation rule: shifts +-pi/2, coefficients +-1/2."""
    return rule_for_frequencies((1.0,))

four_term_rule

four_term_rule() -> ShiftRule

Controlled rotations: generator eigenvalues {0, 0, +-1/2} => frequencies {1/2, 1}.

Source code in src/qmlkit/gradients/rules.py

def four_term_rule() -> ShiftRule:
    """Controlled rotations: generator eigenvalues {0, 0, +-1/2} => frequencies {1/2, 1}."""
    return rule_for_frequencies((0.5, 1.0))

rule_for_gate

rule_for_gate(gate: str) -> ShiftRule

Look the rule up from the gate's declared generator frequencies.

This per-gate lookup is the whole point: a circuit mixing ry with crz needs two different rules, and applying one uniform rule to both returns a plausible, wrong gradient with no error raised.

Source code in src/qmlkit/gradients/rules.py

def rule_for_gate(gate: str) -> ShiftRule:
    """Look the rule up from the gate's declared generator frequencies.

    This per-gate lookup is the whole point: a circuit mixing ``ry`` with ``crz``
    needs two *different* rules, and applying one uniform rule to both returns a
    plausible, wrong gradient with no error raised.
    """
    g = get_gate(gate)
    if not g.is_parametric:
        raise ValueError(f"gate {gate!r} has no parameters to differentiate")
    if not g.is_differentiable:
        raise ValueError(
            f"gate {gate!r} declares no generator frequencies, so no exact shift rule "
            "can be derived. Register it with frequencies=(...) to make it differentiable."
        )
    return rule_for_frequencies(g.frequencies)

second_derivative_rule

second_derivative_rule(
    frequencies: tuple[float, ...] = (1.0,),
) -> ShiftRule

Diagonal Hessian rule. For one frequency: f'' = (f(theta+pi) - f(theta)) / 2.

Source code in src/qmlkit/gradients/rules.py

def second_derivative_rule(frequencies: tuple[float, ...] = (1.0,)) -> ShiftRule:
    """Diagonal Hessian rule. For one frequency: f'' = (f(theta+pi) - f(theta)) / 2."""
    if frequencies != (1.0,):
        raise NotImplementedError(
            "the closed-form second-derivative rule is implemented for single-frequency "
            "(Pauli) generators only; apply the first-order rule twice otherwise"
        )
    return ShiftRule(shifts=(np.pi,), coeffs=(0.5,), unshifted_coeff=-0.5)

qmlkit.gradients.parameter_shift

parameter_shift

The parameter-shift rule — exact gradients from measurements alone.

Two things here are easy to get wrong, and both fail silently — a plausible number, no exception:

Per-gate rules. The shift rule is a property of the gate's generator, not of the call. A circuit mixing ry (one frequency, two-term rule) with crz (two frequencies, four-term rule) needs both. We look the rule up per slot.

Per-occurrence shifting. When one logical parameter drives several gates — weight tying, as in a QCNN's shared convolution block — the derivative is the sum over occurrences, each shifted on its own. Shifting them together computes something else. The slot abstraction in :mod:qmlkit.core.ir makes this fall out naturally: several slots simply map back to the same parameter index.

param_shift_grad

param_shift_grad(
    f_slots: SlotFn,
    spec: CircuitSpec,
    theta: Sequence[float],
    rules: dict[int, ShiftRule] | None = None,
    f0: float | None = None,
) -> NDArray[Any]

Exact gradient of f with respect to the logical parameter vector.

Parameters:

Name	Type	Description	Default
f_slots	SlotFn	Evaluates the circuit given slot angles — one angle per parameterised gate site. Taking slot angles rather than the logical vector is what makes per-occurrence shifting expressible at all.	required
spec	CircuitSpec	Supplies the slot map and each slot's gate (hence its shift rule).	required
theta	Sequence[float]	The logical parameter vector.	required
rules	dict[int, ShiftRule] | None	Optional {slot_index: ShiftRule} override.	None
f0	float | None	The unshifted value, if you already have it. Unused by the standard two-term rule; required by rules with needs_unshifted.	None

Source code in src/qmlkit/gradients/parameter_shift.py

def param_shift_grad(
    f_slots: SlotFn,
    spec: CircuitSpec,
    theta: Sequence[float],
    rules: dict[int, ShiftRule] | None = None,
    f0: float | None = None,
) -> npt.NDArray[Any]:
    """Exact gradient of ``f`` with respect to the logical parameter vector.

    Parameters
    ----------
    f_slots
        Evaluates the circuit given **slot angles** — one angle per parameterised
        gate site. Taking slot angles rather than the logical vector is what makes
        per-occurrence shifting expressible at all.
    spec
        Supplies the slot map and each slot's gate (hence its shift rule).
    theta
        The logical parameter vector.
    rules
        Optional ``{slot_index: ShiftRule}`` override.
    f0
        The unshifted value, if you already have it. Unused by the standard
        two-term rule; required by rules with ``needs_unshifted``.
    """
    angles = spec.bind_slots(theta)
    slots = spec.slots()
    grad = np.zeros(spec.n_params, dtype=float)

    for i, slot in enumerate(slots):
        rule = (rules or {}).get(i) or rule_for_gate(slot.gate)
        total = 0.0
        for shift, coeff in zip(rule.shifts, rule.coeffs, strict=False):
            shifted = angles.copy()
            shifted[i] += shift
            total += coeff * f_slots(shifted)
        if rule.needs_unshifted:
            base = f0 if f0 is not None else f_slots(angles)
            total += rule.unshifted_coeff * base
        # Chain rule for a rescaled reference, and — crucially — ``+=``, which is
        # what sums a weight-tied parameter's several occurrences.
        grad[slot.ref.index] += total * slot.ref.scale

    return grad

param_shift_grad_circuit

param_shift_grad_circuit(
    spec: CircuitSpec,
    theta: Sequence[float],
    obs: Observable | None = None,
    shots: int | None = None,
    backend: BackendLike = None,
    seed: int | None = None,
) -> NDArray[Any]

Convenience wrapper: parameter-shift gradient of <obs> for a circuit.

Source code in src/qmlkit/gradients/parameter_shift.py

def param_shift_grad_circuit(
    spec: CircuitSpec,
    theta: Sequence[float],
    obs: Observable | None = None,
    shots: int | None = None,
    backend: BackendLike = None,
    seed: int | None = None,
) -> npt.NDArray[Any]:
    """Convenience wrapper: parameter-shift gradient of ``<obs>`` for a circuit."""
    obs = Z(0) if obs is None else obs
    be = get_backend(backend)

    def f_slots(angles: npt.NDArray[Any]) -> float:
        return be.expectation(spec.with_slot_angles(angles), obs, shots, seed)

    return param_shift_grad(f_slots, spec, theta)

grad_circuit_cost

grad_circuit_cost(spec: CircuitSpec) -> int

Circuit evaluations for one full parameter-shift gradient.

Sums each slot's real rule cost. This is deliberately not a flat 2P: a controlled rotation costs four evaluations, not two, and a weight-tied parameter costs its rule once per occurrence.

Source code in src/qmlkit/gradients/parameter_shift.py

def grad_circuit_cost(spec: CircuitSpec) -> int:
    """Circuit evaluations for one full parameter-shift gradient.

    Sums each slot's *real* rule cost. This is deliberately not a flat ``2P``:
    a controlled rotation costs four evaluations, not two, and a weight-tied
    parameter costs its rule once per occurrence.
    """
    return sum(rule_for_gate(s.gate).n_evaluations for s in spec.slots())

finite_diff_grad

finite_diff_grad(
    f: Callable[[NDArray[Any]], float],
    theta: Sequence[float],
    eps: float = 1e-06,
    mode: str = "central",
) -> NDArray[Any]

Finite differences. For debugging and tests only — never for training.

Carries an O(eps**2) bias and amplifies any sampling noise by 1/eps, which is why there is no good choice of eps on a noisy device.

Source code in src/qmlkit/gradients/parameter_shift.py

def finite_diff_grad(
    f: Callable[[npt.NDArray[Any]], float],
    theta: Sequence[float],
    eps: float = 1e-6,
    mode: str = "central",
) -> npt.NDArray[Any]:
    """Finite differences. For debugging and tests only — never for training.

    Carries an ``O(eps**2)`` bias and amplifies any sampling noise by ``1/eps``,
    which is why there is no good choice of ``eps`` on a noisy device.
    """
    arr = np.asarray(theta, dtype=float).ravel()
    grad = np.zeros_like(arr)
    base = f(arr) if mode == "forward" else None
    for k in range(arr.size):
        plus = arr.copy()
        plus[k] += eps
        if mode == "central":
            minus = arr.copy()
            minus[k] -= eps
            grad[k] = (f(plus) - f(minus)) / (2 * eps)
        elif mode == "forward":
            assert base is not None
            grad[k] = (f(plus) - base) / eps
        else:
            raise unknown("mode", mode, ("central", "forward"))
    return grad

qmlkit.gradients.adjoint

adjoint

Adjoint differentiation — exact gradients in one backward pass.

Parameter-shift needs 2P circuit evaluations because a real device can only be run, never inspected. On a simulator the state is right there, so the whole gradient comes out of a single forward and a single backward sweep — independent of P.

The sweep. With :math:|\psi_j\rangle = U_j \cdots U_1 |0\rangle and :math:E = \langle\psi_n| O |\psi_n\rangle,

.. math:: \frac{\partial E}{\partial\theta_k} = 2\,\mathrm{Re}\,\langle\lambda_k|\,\partial_k U_k\,|\psi_{k-1}\rangle, \qquad \lambda_k = (U_n\cdots U_{k+1})^\dagger O |\psi_n\rangle

so walking backwards and undoing one gate at a time keeps both states current at O(1) extra memory.

This is exact — the gate derivatives are closed-form, not finite differences — but it is simulator-only, because no device will hand you an amplitude. Since 0.x is simulator-only anyway, it is the right default for training; parameter-shift remains what the course teaches, what validates this, and what keeps the library hardware-ready.

supports_adjoint

supports_adjoint(
    spec: CircuitSpec, backend: Backend | str | None = None
) -> bool

True if every parameterised gate in spec has a closed-form derivative.

Source code in src/qmlkit/gradients/adjoint.py

def supports_adjoint(spec: CircuitSpec, backend: Backend | str | None = None) -> bool:
    """True if every parameterised gate in ``spec`` has a closed-form derivative."""
    from qmlkit.core.gates import get_gate

    be = get_backend(backend)
    if not be.supports_statevector:
        return False
    return all(get_gate(s.gate).has_derivative for s in spec.slots())

adjoint_grad

adjoint_grad(
    spec: CircuitSpec,
    theta: NDArray[Any],
    obs: Observable | None = None,
    backend: Backend | str | None = None,
) -> NDArray[Any]

Exact gradient of <obs> with respect to the logical parameter vector.

One forward pass and one backward pass, whatever P is. Weight-tied parameters accumulate across their occurrences, exactly as parameter-shift does.

Source code in src/qmlkit/gradients/adjoint.py

def adjoint_grad(
    spec: CircuitSpec,
    theta: npt.NDArray[Any],
    obs: Observable | None = None,
    backend: Backend | str | None = None,
) -> npt.NDArray[Any]:
    """Exact gradient of ``<obs>`` with respect to the logical parameter vector.

    One forward pass and one backward pass, whatever ``P`` is. Weight-tied
    parameters accumulate across their occurrences, exactly as parameter-shift does.
    """
    obs = Z(0) if obs is None else obs
    be = get_backend(backend)
    if not be.supports_statevector:
        raise ValueError(
            f"the {be.name!r} backend has no statevector, so it cannot differentiate by "
            'the adjoint method; use grad_method="parameter-shift"'
        )

    theta = np.asarray(theta, dtype=float).ravel()
    slots = spec.slots()
    slot_angles = spec.bind_slots(theta)
    n = spec.n_qubits
    shape = (2,) * n

    # forward: run the bound circuit once
    psi = be.statevector(spec.with_slot_angles(slot_angles)).reshape(shape)
    lam = _apply_observable(psi, obs, n)

    grad = np.zeros(spec.n_params, dtype=float)
    slot_of_op = {s.op_index: i for i, s in enumerate(slots)}

    # backward: undo one gate at a time, reading off each parameter's contribution
    cursor = len(slot_angles)
    for op_index in range(len(spec.ops) - 1, -1, -1):
        op = spec.ops[op_index]
        slot_i = slot_of_op.get(op_index)
        if slot_i is not None:
            cursor -= 1
            angles = [float(slot_angles[cursor])]
        else:
            angles = [float(p) for p in op.params if not isinstance(p, ParamRef)]

        u = gate_matrix(op.gate, tuple(angles))
        psi = _apply(psi, u.conj().T, op.qubits)  # psi is now the state before this gate

        if slot_i is not None:
            du = gate_derivative(op.gate, tuple(angles))
            mu = _apply(psi, du, op.qubits)
            contribution = 2.0 * float(np.real(np.vdot(lam, mu)))
            ref = slots[slot_i].ref
            grad[ref.index] += contribution * ref.scale  # += ties occurrences together

        lam = _apply(lam, u.conj().T, op.qubits)

    return grad

adjoint_grad_terms

adjoint_grad_terms(
    spec: CircuitSpec,
    theta: NDArray[Any],
    obs: Observable,
    backend: Backend | str | None = None,
) -> dict[PauliString, NDArray[Any]]

Per-term gradients, for diagnosing which observable term drives a parameter.

Source code in src/qmlkit/gradients/adjoint.py

def adjoint_grad_terms(
    spec: CircuitSpec,
    theta: npt.NDArray[Any],
    obs: Observable,
    backend: Backend | str | None = None,
) -> dict[PauliString, npt.NDArray[Any]]:  # pragma: no cover - convenience for diagnostics
    """Per-term gradients, for diagnosing which observable term drives a parameter."""
    return {t: adjoint_grad(spec, theta, t, backend) for t in as_sum(obs).terms}

qmlkit.gradients.hadamard

hadamard

Hadamard-test gradient — one circuit per parameter instead of two.

For :math:U_k = e^{-i\theta_k P_k/2} inserted at position :math:k, write :math:|\varphi\rangle for the state you get by applying :math:P_k right after :math:U_k. Then

.. math:: \partial_k E = -\,\mathrm{Im}\,\langle\varphi| O |\psi\rangle

and that imaginary part is exactly what a Hadamard test reads out. Put an ancilla in :math:|+\rangle, run the circuit with a controlled :math:P_k inserted after gate :math:k, and measure :math:\langle Y_a \otimes O\rangle:

.. code-block:: text

|+>_a ----------●----------  <Y_a (x) O>
|0>_n -- U_1..U_k -- P_k -- U_{k+1}..U_L --

Half the circuits of parameter-shift, and unlike adjoint it is a real measurement, so it is valid on hardware. The trade is an ancilla that must couple to every wire the generator touches, plus controlled gates. On real devices that routing cost usually outweighs the saved circuits, which is why parameter-shift stays the default there.

On a simulator the halving does show up in wall-clock: 5 qubits, P=120, two-term observable gives 404 ms against parameter-shift's 823 ms. Getting that required passing the lifted observable as a single sum rather than term by term — the ancilla already doubles the statevector, so re-preparing it per term would have handed the whole advantage straight back.

Only single-qubit Pauli rotations are supported: their generators are Paulis, so the controlled form is a cx/cy/cz. A controlled rotation's generator is not a Pauli, and this method refuses it rather than guessing.

supports_hadamard_grad

supports_hadamard_grad(spec: CircuitSpec) -> bool

True if every parameterised gate has a plain Pauli generator.

Source code in src/qmlkit/gradients/hadamard.py

def supports_hadamard_grad(spec: CircuitSpec) -> bool:
    """True if every parameterised gate has a plain Pauli generator."""
    return all(s.gate in _GENERATOR for s in spec.slots())

hadamard_grad

hadamard_grad(
    spec: CircuitSpec,
    theta: NDArray[Any],
    obs: Observable | None = None,
    backend: BackendLike = None,
    shots: int | None = None,
    seed: int | None = None,
) -> NDArray[Any]

Exact gradient using one extra qubit and one circuit per parameter.

Source code in src/qmlkit/gradients/hadamard.py

def hadamard_grad(
    spec: CircuitSpec,
    theta: npt.NDArray[Any],
    obs: Observable | None = None,
    backend: BackendLike = None,
    shots: int | None = None,
    seed: int | None = None,
) -> npt.NDArray[Any]:
    """Exact gradient using one extra qubit and one circuit per parameter."""
    obs = Z(0) if obs is None else obs
    if not supports_hadamard_grad(spec):
        offenders = sorted({s.gate for s in spec.slots() if s.gate not in _GENERATOR})
        raise ValueError(
            f"the Hadamard-test gradient needs Pauli-generated rotations; {offenders} "
            'are not. Use method="parameter-shift", which handles any declared spectrum.'
        )

    arr = np.asarray(theta, dtype=float).ravel()
    slots = spec.slots()
    slot_angles = spec.bind_slots(arr)
    bound = spec.with_slot_angles(slot_angles)
    n = spec.n_qubits
    ancilla = n

    grad = np.zeros(spec.n_params, dtype=float)
    for slot in slots:
        pauli = _GENERATOR[slot.gate]
        target = spec.ops[slot.op_index].qubits[0]

        ops: list[Op] = [Op("h", (ancilla,))]
        for j, op in enumerate(bound.ops):
            ops.append(op)
            if j == slot.op_index:
                # controlled generator, inserted immediately after the gate
                ops.append(Op(_CONTROLLED[pauli], (ancilla, target)))
        probe = CircuitSpec(n + 1, tuple(ops), 0)

        # <Y_ancilla (x) O> reads off the imaginary part we need. Pass the whole sum
        # in one call: the backend prepares the state once and accumulates every term
        # from it, so a k-term observable still costs a single circuit, not k of them.
        lifted = PauliSum(
            tuple(PauliString((*t.paulis, (ancilla, "Y")), t.coeff) for t in as_sum(obs).terms)
        )
        value = expval(probe, lifted, shots=shots, backend=backend, seed=seed)
        grad[slot.ref.index] += value * slot.ref.scale  # += sums tied occurrences

    return grad

hadamard_grad_cost

hadamard_grad_cost(spec: CircuitSpec) -> int

Circuits per gradient: one per parameterised slot.

Source code in src/qmlkit/gradients/hadamard.py

def hadamard_grad_cost(spec: CircuitSpec) -> int:
    """Circuits per gradient: one per parameterised slot."""
    return len(spec.slots())

qmlkit.gradients.spsa

spsa

SPSA — a gradient estimate in two evaluations, whatever P is.

Parameter-shift costs 2P circuits. SPSA perturbs every parameter at once along a random :math:\pm 1 direction and costs two, forever:

.. math:: \hat{g} = \frac{f(\theta + c\Delta) - f(\theta - c\Delta)}{2c}\,\Delta^{-1}

The estimate is noisy but unbiased in expectation, and stochastic optimisers tolerate that well — which is why it is the standard answer once P gets large enough that 2P circuits per step stops being affordable.

The decay schedules are Spall's. A, the stability constant, is the one people leave out: without it the first steps are far too large, and the run diverges before the schedule has a chance to settle it. Roughly 10 % of the planned iteration count is the usual choice.

SPSASchedule

SPSASchedule(
    a: float = 0.2,
    c: float = 0.1,
    A: float | None = None,
    alpha: float = 0.602,
    gamma: float = 0.101,
    n_iterations: int = 100,
)

Spall's decay schedules for the step size and the perturbation size.

Source code in src/qmlkit/gradients/spsa.py

def __init__(
    self,
    a: float = 0.2,
    c: float = 0.1,
    A: float | None = None,
    alpha: float = 0.602,
    gamma: float = 0.101,
    n_iterations: int = 100,
) -> None:
    self.a = a
    self.c = c
    # the stability constant most write-ups omit; without it the early steps
    # are the largest, which is backwards
    self.A = A if A is not None else max(1.0, 0.1 * n_iterations)
    self.alpha = alpha
    self.gamma = gamma

spsa_grad

spsa_grad(
    f: LossFn,
    theta: ArrayLike,
    c: float = 0.1,
    n_avg: int = 1,
    seed: int | None = None,
    rng: Generator | None = None,
) -> NDArray[Any]

A stochastic gradient estimate from 2 * n_avg evaluations.

n_avg averages several random directions, trading evaluations for variance — still constant in P.

Source code in src/qmlkit/gradients/spsa.py

def spsa_grad(
    f: LossFn,
    theta: ArrayLike,
    c: float = 0.1,
    n_avg: int = 1,
    seed: int | None = None,
    rng: np.random.Generator | None = None,
) -> npt.NDArray[Any]:
    """A stochastic gradient estimate from ``2 * n_avg`` evaluations.

    ``n_avg`` averages several random directions, trading evaluations for variance
    — still constant in ``P``.
    """
    arr = np.asarray(theta, dtype=float).ravel()
    if c <= 0:
        raise ValueError("c must be positive")
    if n_avg < 1:
        raise ValueError("n_avg must be at least 1")
    generator = rng if rng is not None else np.random.default_rng(seed)

    total = np.zeros_like(arr)
    for _ in range(n_avg):
        delta = generator.choice([-1.0, 1.0], size=arr.shape)
        plus = float(f(arr + c * delta))
        minus = float(f(arr - c * delta))
        total += (plus - minus) / (2.0 * c) * delta  # delta^-1 == delta for +-1
    return total / n_avg

spsa_step

spsa_step(
    f: LossFn,
    theta: ArrayLike,
    k: int,
    schedule: SPSASchedule | None = None,
    rng: Generator | None = None,
) -> NDArray[Any]

One SPSA update at iteration k.

Source code in src/qmlkit/gradients/spsa.py

def spsa_step(
    f: LossFn,
    theta: ArrayLike,
    k: int,
    schedule: SPSASchedule | None = None,
    rng: np.random.Generator | None = None,
) -> npt.NDArray[Any]:
    """One SPSA update at iteration ``k``."""
    sched = schedule or SPSASchedule()
    arr = np.asarray(theta, dtype=float).ravel()
    g = spsa_grad(f, arr, c=sched.perturbation(k), rng=rng)
    return arr - sched.step_size(k) * g

minimize_spsa

minimize_spsa(
    f: LossFn,
    theta0: ArrayLike,
    n_iterations: int = 100,
    schedule: SPSASchedule | None = None,
    seed: int | None = None,
    callback: Callable[[int, NDArray[Any], float], None]
    | None = None,
) -> tuple[NDArray[Any], list[float]]

Minimise f with SPSA. Returns the final parameters and the loss history.

Two evaluations per iteration regardless of how many parameters there are.

Source code in src/qmlkit/gradients/spsa.py

def minimize_spsa(
    f: LossFn,
    theta0: ArrayLike,
    n_iterations: int = 100,
    schedule: SPSASchedule | None = None,
    seed: int | None = None,
    callback: Callable[[int, npt.NDArray[Any], float], None] | None = None,
) -> tuple[npt.NDArray[Any], list[float]]:
    """Minimise ``f`` with SPSA. Returns the final parameters and the loss history.

    Two evaluations per iteration regardless of how many parameters there are.
    """
    sched = schedule or SPSASchedule(n_iterations=n_iterations)
    rng = np.random.default_rng(seed)
    # Annotated explicitly: newer NumPy stubs give `.ravel()` a *shape-typed* 1-D
    # result, and reassigning a general-shape array to it then fails strict type-checking.
    # Only CI's Python 3.10 NumPy tripped on this; two other NumPy generations did not.
    theta: npt.NDArray[Any] = np.asarray(theta0, dtype=float).ravel().copy()
    history: list[float] = []
    for k in range(n_iterations):
        value = float(f(theta))
        history.append(value)
        if callback is not None:
            callback(k, theta, value)
        theta = spsa_step(f, theta, k, sched, rng)
    history.append(float(f(theta)))
    return theta, history