Analysis

Measuring an ansatz rather than asserting things about it: expressibility, entanglement, spectra, geometry.

qmlkit.diagnostics

diagnostics

Whether the model you just built is quietly broken.

In most libraries a mistake raises. In quantum machine learning it returns a number — right shape, right range, entirely plausible, and wrong:

• A re-uploading model whose trainable block commutes with its encoding trains, converges, and reaches one Fourier frequency instead of the L it was designed for. Every weight in it is a phase shift.
• An ansatz with a parameter the circuit cannot feel fits exactly as well without it. The optimiser reports no difficulty, because there is none.
• A kernel that has concentrated gives every pair of inputs the same similarity, and still produces a Gram matrix, an SVM, and an accuracy.

None of these are exceptions to catch. They are properties to measure, and this module measures them.

The checks are deliberately decisive rather than exhaustive. Each one has a threshold that separates "wrong" from "unusual" with room to spare, because a diagnostic that cries wolf is one nobody runs twice. Where the check is exact — a parameter that cannot change the state is dead, full stop — the threshold is machine epsilon. Where it is statistical, the finding says what was measured, so the number can be argued with.

>>> import qmlkit as qk
>>> report = qk.diagnose(qk.hardware_efficient(3, 2))
>>> bool(report)
False

A report is falsy when it found nothing, so if qk.diagnose(model): ... reads the way it should. Findings carry a stable code to branch on, a message saying what is wrong, and a fix naming the edit that resolves it.

Finding dataclass

Finding(
    code: str,
    severity: str,
    message: str,
    fix: str = "",
    value: float | None = None,
)

One thing that is wrong, why it matters, and the edit that fixes it.

Report dataclass

Report(subject: str, findings: tuple[Finding, ...] = ())

Everything :func:diagnose found, worst first.

Falsy when empty, so it can be tested directly. Iterating yields :class:Finding objects; codes is the flat list to assert against.

diagnose

diagnose(
    subject: object,
    X: Any = None,
    y: Any = None,
    *,
    obs: Observable | None = None,
    n_samples: int = 30,
    probes: int = 3,
    seed: int | None = 0,
    backend: BackendLike = None,
    shots: int | None = None,
    n_qubits: int | None = None,
) -> Report

Check a model or a Gram matrix for the failures that do not raise.

Parameters:

Name	Type	Description	Default
subject	object	An :class:~qmlkit.ansatz.library.Ansatz, anything holding one (a QuantumLayer, VQC, VQRegressor, or an nn.Sequential containing one), or a square Gram matrix.	required
X	Any	Optional. Given both, and a trained model holding a QuantumLayer, the structural checks are joined by one that needs data: whether the quantum layer earned its place, or the classical layers around it are carrying the model. Classification targets only -- the probe is a classifier.	None
y	Any	Optional. Given both, and a trained model holding a QuantumLayer, the structural checks are joined by one that needs data: whether the quantum layer earned its place, or the classical layers around it are carrying the model. Classification targets only -- the probe is a classifier.	None
obs	Observable | None	Observable for the trainability probe. Defaults to Z(0), matching :func:~qmlkit.metrics.barren_plateau_scan.	None
n_samples	int	Sample count for the statistical checks — entanglement and gradient variance. The exact checks ignore it.	30
probes	int	Random points at which to test whether a parameter can move the state. A parameter is dead if it moves nothing at any of them; three is already conclusive, since the points are random and the test is exact.	3
shots	int | None	Gram matrices only. shots enables the check for whether the signal survives sampling noise; n_qubits enables the comparison against the 2^-n concentration scale.	None
n_qubits	int | None	Gram matrices only. shots enables the check for whether the signal survives sampling noise; n_qubits enables the comparison against the 2^-n concentration scale.	None

Returns:

Type	Description
Report	Falsy when nothing was found. Sorted worst first.

Examples:

>>> import qmlkit as qk
>>> healthy = qk.diagnose(qk.hardware_efficient(3, 2))
>>> bool(healthy)
False

A model whose weights share the encoding's generator is the trap the re-uploading literature warns about, and it is silent without this:

>>> from qmlkit.ansatz import Ansatz, EncodingLayer, RotationLayer, repeat
>>> fmap = qk.AngleFeatureMap(1, rotation="ry")
>>> block = EncodingLayer(fmap) + RotationLayer("ry")
>>> broken = qk.diagnose(Ansatz(1, repeat(3, block), n_inputs=1))
>>> "ENCODING_COMMUTES" in broken.codes
True

Source code in src/qmlkit/diagnostics.py

def diagnose(
    subject: object,
    X: Any = None,
    y: Any = None,
    *,
    obs: Observable | None = None,
    n_samples: int = 30,
    probes: int = 3,
    seed: int | None = 0,
    backend: BackendLike = None,
    shots: int | None = None,
    n_qubits: int | None = None,
) -> Report:
    """Check a model or a Gram matrix for the failures that do not raise.

    Parameters
    ----------
    subject
        An :class:`~qmlkit.ansatz.library.Ansatz`, anything holding one (a
        ``QuantumLayer``, ``VQC``, ``VQRegressor``, or an ``nn.Sequential``
        containing one), or a square Gram matrix.
    X, y
        Optional. Given both, and a *trained* model holding a ``QuantumLayer``, the
        structural checks are joined by one that needs data: whether the quantum
        layer earned its place, or the classical layers around it are carrying the
        model. Classification targets only -- the probe is a classifier.
    obs
        Observable for the trainability probe. Defaults to ``Z(0)``, matching
        :func:`~qmlkit.metrics.barren_plateau_scan`.
    n_samples
        Sample count for the statistical checks — entanglement and gradient
        variance. The exact checks ignore it.
    probes
        Random points at which to test whether a parameter can move the state. A
        parameter is dead if it moves nothing at any of them; three is already
        conclusive, since the points are random and the test is exact.
    shots, n_qubits
        Gram matrices only. ``shots`` enables the check for whether the signal
        survives sampling noise; ``n_qubits`` enables the comparison against the
        ``2^-n`` concentration scale.

    Returns
    -------
    Report
        Falsy when nothing was found. Sorted worst first.

    Examples
    --------
    >>> import qmlkit as qk
    >>> healthy = qk.diagnose(qk.hardware_efficient(3, 2))
    >>> bool(healthy)
    False

    A model whose weights share the encoding's generator is the trap the
    re-uploading literature warns about, and it is silent without this:

    >>> from qmlkit.ansatz import Ansatz, EncodingLayer, RotationLayer, repeat
    >>> fmap = qk.AngleFeatureMap(1, rotation="ry")
    >>> block = EncodingLayer(fmap) + RotationLayer("ry")
    >>> broken = qk.diagnose(Ansatz(1, repeat(3, block), n_inputs=1))
    >>> "ENCODING_COMMUTES" in broken.codes
    True
    """
    if isinstance(subject, np.ndarray) or (
        isinstance(subject, list | tuple) and subject and isinstance(subject[0], list | tuple)
    ):
        gram = np.asarray(subject, dtype=float)
        return Report(
            f"Gram matrix {gram.shape[0]}x{gram.shape[0]}",
            _sorted(_diagnose_kernel(gram, n_qubits=n_qubits, shots=shots)),
        )

    ansatz = _find_ansatz(subject)
    if ansatz is None:
        raise TypeError(
            f"diagnose() cannot inspect a {type(subject).__name__}. It takes an Ansatz, "
            "anything holding one (QuantumLayer, VQC, VQRegressor, or an nn.Sequential "
            "containing one), or a Gram matrix as a square array."
        )
    feature_map = _find_feature_map(subject)
    prefix = feature_map.build_parametric(offset=0) if feature_map is not None else None
    findings = _diagnose_ansatz(
        ansatz,
        obs=obs,
        n_samples=n_samples,
        probes=probes,
        seed=seed,
        backend=backend,
        prefix=prefix,
        observables=_find_observables(subject),
    )
    if X is not None and y is not None:
        findings += _diagnose_contribution(subject, X, y, seed=seed)
    label = f"{type(subject).__name__} ({ansatz.name})" if subject is not ansatz else ansatz.name
    subject_line = f"{label} on {ansatz.n_qubits} qubits"
    _, substituted = _structural_backend(backend)
    if substituted is not None:
        subject_line += (
            f" [structure checked on the numpy reference: {substituted!r} has no statevector]"
        )
    return Report(subject_line, _sorted(findings))

qmlkit.metrics

metrics

Does this ansatz stand a chance? — expressibility, entanglement, trainability.

Choosing an ansatz by eye is guesswork. These are the four numbers the literature actually uses, and they pull against each other:

• Expressibility — how close the ansatz's state distribution gets to Haar-random. Measured as KL(fidelities || Haar); smaller is more expressive.
• Entangling capability — the mean Meyer–Wallach Q; 0 is a product state, 1 is maximally entangled.
• Trainability — the variance of the gradient. It collapses roughly as 2^-n for deep circuits, which is the barren plateau: gradients vanish faster than any shot budget can resolve them.
• Generalization — how much data the model needs, growing with trainable gates.

More expressibility costs trainability. That trade is the whole design problem, and :class:AnsatzReport puts both numbers side by side so the choice is informed.

AnsatzReport dataclass

AnsatzReport(
    ansatz: Ansatz,
    n_samples: int = 300,
    seed: int | None = 0,
    backend: BackendLike = None,
    results: dict[str, object] = dict(),
)

Expressibility, entanglement, depth, cost and trainability in one call.

print(AnsatzReport(qk.hardware_efficient(4, 2)))

haar_fidelity_pdf

haar_fidelity_pdf(
    f: NDArray[Any], n_qubits: int
) -> NDArray[Any]

Haar-random fidelity density: (N-1)(1-F)^(N-2) for N = 2^n.

Source code in src/qmlkit/metrics.py

def haar_fidelity_pdf(f: npt.NDArray[Any], n_qubits: int) -> npt.NDArray[Any]:
    """Haar-random fidelity density: ``(N-1)(1-F)^(N-2)`` for ``N = 2^n``."""
    n_states = 2**n_qubits
    fid = np.clip(np.asarray(f, dtype=float), 0.0, 1.0)
    return (n_states - 1) * (1.0 - fid) ** (n_states - 2)

fidelity_samples

fidelity_samples(
    ansatz: Ansatz,
    n_samples: int = 2000,
    seed: int | None = None,
    backend: BackendLike = None,
) -> NDArray[Any]

Fidelities between pairs of states from independently sampled parameters.

Source code in src/qmlkit/metrics.py

def fidelity_samples(
    ansatz: Ansatz,
    n_samples: int = 2000,
    seed: int | None = None,
    backend: BackendLike = None,
) -> npt.NDArray[Any]:
    """Fidelities between pairs of states from independently sampled parameters."""
    require_statevector(backend, "fidelity_samples")
    rng = np.random.default_rng(seed)
    out = np.empty(n_samples, dtype=float)
    for i in range(n_samples):
        a = statevector(ansatz.build(rng.uniform(-np.pi, np.pi, ansatz.n_params)), backend=backend)
        b = statevector(ansatz.build(rng.uniform(-np.pi, np.pi, ansatz.n_params)), backend=backend)
        out[i] = abs(np.vdot(a, b)) ** 2
    return out

expressibility

expressibility(
    ansatz: Ansatz,
    n_samples: int = 2000,
    n_bins: int = 75,
    seed: int | None = None,
    backend: BackendLike = None,
) -> float

KL(ansatz fidelities || Haar). Smaller is more expressive; 0 is Haar.

Note the direction — it is a divergence from Haar, so a low number means the ansatz reaches as much of state space as a random circuit would.

Source code in src/qmlkit/metrics.py

def expressibility(
    ansatz: Ansatz,
    n_samples: int = 2000,
    n_bins: int = 75,
    seed: int | None = None,
    backend: BackendLike = None,
) -> float:
    """``KL(ansatz fidelities || Haar)``. **Smaller is more expressive**; 0 is Haar.

    Note the direction — it is a divergence *from* Haar, so a low number means the
    ansatz reaches as much of state space as a random circuit would.
    """
    # named here as well as in fidelity_samples, so the refusal names the function
    # the caller actually invoked
    require_statevector(backend, "expressibility")
    fids = fidelity_samples(ansatz, n_samples, seed, backend)
    edges = np.linspace(0.0, 1.0, n_bins + 1)
    observed, _ = np.histogram(fids, bins=edges)
    p = observed / observed.sum()

    centres = (edges[:-1] + edges[1:]) / 2
    haar = haar_fidelity_pdf(centres, ansatz.n_qubits)
    q = haar / haar.sum()

    mask = p > 0
    return float(np.sum(p[mask] * np.log(p[mask] / np.clip(q[mask], 1e-300, None))))

meyer_wallach

meyer_wallach(
    state: NDArray[Any], n_qubits: int | None = None
) -> float

Meyer–Wallach Q = 2(1 - (1/n) sum_k Tr rho_k^2).

0 for any product state, 1 for a maximally entangled one.

Source code in src/qmlkit/metrics.py

def meyer_wallach(state: npt.NDArray[Any], n_qubits: int | None = None) -> float:
    """Meyer–Wallach ``Q = 2(1 - (1/n) sum_k Tr rho_k^2)``.

    0 for any product state, 1 for a maximally entangled one.
    """
    psi = np.asarray(state, dtype=complex).ravel()
    n = n_qubits if n_qubits is not None else int(np.log2(psi.size))
    return float(2.0 * (1.0 - np.mean([purity(psi, [q], n) for q in range(n)])))

entangling_capability

entangling_capability(
    ansatz: Ansatz,
    n_samples: int = 200,
    seed: int | None = None,
    backend: BackendLike = None,
) -> float

Mean Meyer–Wallach Q over randomly sampled parameters.

Source code in src/qmlkit/metrics.py

def entangling_capability(
    ansatz: Ansatz,
    n_samples: int = 200,
    seed: int | None = None,
    backend: BackendLike = None,
) -> float:
    """Mean Meyer–Wallach ``Q`` over randomly sampled parameters."""
    require_statevector(backend, "entangling_capability")
    rng = np.random.default_rng(seed)
    vals = [
        meyer_wallach(
            statevector(ansatz.build(rng.uniform(-np.pi, np.pi, ansatz.n_params)), backend=backend),
            ansatz.n_qubits,
        )
        for _ in range(n_samples)
    ]
    return float(np.mean(vals))

gradient_variance

gradient_variance(
    ansatz: Ansatz,
    obs: Observable | None = None,
    n_samples: int = 100,
    param_index: int = 0,
    seed: int | None = None,
    backend: BackendLike = None,
) -> float

Variance of one parameter's gradient over random initialisations.

This is the barren-plateau probe: if it falls exponentially with width, no realistic shot budget will resolve the gradient.

It probes one parameter, and which one matters. The default param_index=0 is a leading rotation on several stock ansaetze, and a leading Rz on |0> has a gradient of exactly zero against Z - not because the ansatz is untrainable but because that parameter does nothing. A variance at machine zero is reported with a warning saying so, because the number alone is indistinguishable from a genuine plateau. :func:~qmlkit.diagnostics.diagnose finds dead parameters directly and names their indices.

Source code in src/qmlkit/metrics.py

def gradient_variance(
    ansatz: Ansatz,
    obs: Observable | None = None,
    n_samples: int = 100,
    param_index: int = 0,
    seed: int | None = None,
    backend: BackendLike = None,
) -> float:
    """Variance of one parameter's gradient over random initialisations.

    This is the barren-plateau probe: if it falls exponentially with width, no
    realistic shot budget will resolve the gradient.

    **It probes one parameter, and which one matters.** The default ``param_index=0``
    is a leading rotation on several stock ansaetze, and a leading ``Rz`` on ``|0>``
    has a gradient of exactly zero against ``Z`` - not because the ansatz is
    untrainable but because that parameter does nothing. A variance at machine zero
    is reported with a warning saying so, because the number alone is
    indistinguishable from a genuine plateau. :func:`~qmlkit.diagnostics.diagnose`
    finds dead parameters directly and names their indices.
    """
    from qmlkit.gradients.dispatch import grad

    obs = Z(0) if obs is None else obs
    rng = np.random.default_rng(seed)
    spec = ansatz.build()
    vals = [
        float(
            grad(spec, rng.uniform(-np.pi, np.pi, ansatz.n_params), obs, backend=backend)[
                param_index
            ]
        )
        for _ in range(n_samples)
    ]
    variance = float(np.var(vals))
    if variance < _DEAD_GRADIENT:
        warnings.warn(
            f"gradient variance for parameter {param_index} is {variance:.3e}, which is "
            "machine zero rather than a small number: that parameter does not move this "
            "observable at all. This reads like a barren plateau and is not one. Probe a "
            "different param_index, or run qmlkit.diagnose(ansatz), which names dead "
            "parameters directly.",
            stacklevel=2,
        )
    return variance

barren_plateau_scan

barren_plateau_scan(
    ansatz_factory: Callable[[int], Ansatz],
    qubit_range: Sequence[int],
    obs_factory: Callable[[int], Observable] | None = None,
    n_samples: int = 100,
    seed: int | None = None,
    backend: BackendLike = None,
) -> dict[str, Any]

Gradient variance against qubit count.

obs_factory decides the cost locality, which matters at fixed shallow depth: measured on a 2-layer hardware-efficient ansatz from 2 to 6 qubits, a local Z(0) holds its gradient variance flat (decay 0.98 per qubit) while a global Z^n collapses exponentially (0.56). Depth eventually wins regardless — at L = 2n both decay exponentially — so this reports the measurement rather than asserting a rule.

Source code in src/qmlkit/metrics.py

def barren_plateau_scan(
    ansatz_factory: Callable[[int], Ansatz],
    qubit_range: Sequence[int],
    obs_factory: Callable[[int], Observable] | None = None,
    n_samples: int = 100,
    seed: int | None = None,
    backend: BackendLike = None,
) -> dict[str, Any]:
    """Gradient variance against qubit count.

    ``obs_factory`` decides the **cost locality**, which matters at fixed shallow
    depth: measured on a 2-layer hardware-efficient ansatz from 2 to 6 qubits, a
    local ``Z(0)`` holds its gradient variance flat (decay 0.98 per qubit) while a
    global ``Z^n`` collapses exponentially (0.56). Depth eventually wins regardless
    — at ``L = 2n`` both decay exponentially — so this reports the measurement
    rather than asserting a rule.
    """
    obs_factory = obs_factory or (lambda n: Z(0))
    widths, variances = [], []
    for n in qubit_range:
        widths.append(n)
        variances.append(
            gradient_variance(ansatz_factory(n), obs_factory(n), n_samples, 0, seed, backend)
        )
    decay = _decay_rate(widths, variances)
    return {
        "n_qubits": widths,
        "variance": variances,
        "decay_per_qubit": decay,
        "looks_exponential": bool(decay is not None and decay < 0.75),
    }

fisher_information

fisher_information(
    ansatz: Ansatz,
    X: NDArray[Any],
    theta: NDArray[Any],
    obs: Observable | None = None,
    backend: BackendLike = None,
) -> NDArray[Any]

Classical Fisher information of the model output, averaged over inputs.

This is the classical FIM of the output distribution — the object effective dimension is built on. Not to be confused with the quantum Fisher information, which is 4 x the Fubini–Study metric and is what natural gradient uses.

Source code in src/qmlkit/metrics.py

def fisher_information(
    ansatz: Ansatz,
    X: npt.NDArray[Any],
    theta: npt.NDArray[Any],
    obs: Observable | None = None,
    backend: BackendLike = None,
) -> npt.NDArray[Any]:
    """Classical Fisher information of the model output, averaged over inputs.

    This is the *classical* FIM of the output distribution — the object effective
    dimension is built on. Not to be confused with the quantum Fisher information,
    which is ``4 x`` the Fubini–Study metric and is what natural gradient uses.
    """
    from qmlkit.gradients.dispatch import grad

    obs = Z(0) if obs is None else obs
    spec = ansatz.build()
    p = len(theta)
    total = np.zeros((p, p))
    rows = np.atleast_2d(np.asarray(X, dtype=float))
    for _ in rows:
        g = grad(spec, theta, obs, backend=backend)
        total += np.outer(g, g)
    return total / max(len(rows), 1)

effective_dimension

effective_dimension(
    fisher: NDArray[Any],
    n_samples: int = 1000,
    gamma: float = 1.0,
) -> float

Normalised effective dimension of a model, from its Fisher information.

How many parameters are usefully independent, which is generally far fewer than the raw count. Follows the normalised-Fisher construction rather than the "count eigenvalues above a threshold" shortcut, which is a teaching stand-in.

Source code in src/qmlkit/metrics.py

def effective_dimension(
    fisher: npt.NDArray[Any], n_samples: int = 1000, gamma: float = 1.0
) -> float:
    """Normalised effective dimension of a model, from its Fisher information.

    How many parameters are *usefully* independent, which is generally far fewer
    than the raw count. Follows the normalised-Fisher construction rather than the
    "count eigenvalues above a threshold" shortcut, which is a teaching stand-in.
    """
    f = np.asarray(fisher, dtype=float)
    p = f.shape[0]
    trace = np.trace(f)
    if trace <= 0:
        return 0.0
    f_hat = p * f / trace  # normalised so the trace is p
    kappa = gamma * n_samples / (2 * np.pi * np.log(n_samples))
    eig = np.linalg.eigvalsh(f_hat)
    numerator = float(np.sum(np.log1p(np.clip(kappa * eig, 0, None))))
    denominator = 2.0 * np.log(kappa) if kappa > 1 else 1.0
    return float(np.clip(numerator / denominator, 0.0, p))

generalization_bound

generalization_bound(
    n_trainable_gates: int,
    n_samples: int,
    with_log: bool = True,
) -> float

Expected generalization gap, O(sqrt(T log T / N)) (Caro et al. 2022).

with_log=False drops the log factor for the simplified sqrt(T/N) form often quoted — easier to reason about, but not the actual bound.

Source code in src/qmlkit/metrics.py

def generalization_bound(n_trainable_gates: int, n_samples: int, with_log: bool = True) -> float:
    r"""Expected generalization gap, ``O(sqrt(T log T / N))`` (Caro et al. 2022).

    ``with_log=False`` drops the log factor for the simplified ``sqrt(T/N)`` form
    often quoted — easier to reason about, but not the actual bound.
    """
    if n_samples <= 0:
        raise ValueError("n_samples must be positive")
    t = max(int(n_trainable_gates), 1)
    numerator = t * np.log(t) if with_log and t > 1 else t
    return float(np.sqrt(numerator / n_samples))

samples_for_gap

samples_for_gap(
    n_trainable_gates: int,
    gap: float,
    with_log: bool = True,
) -> int

Invert the bound: how many examples to reach a target generalization gap.

Source code in src/qmlkit/metrics.py

def samples_for_gap(n_trainable_gates: int, gap: float, with_log: bool = True) -> int:
    """Invert the bound: how many examples to reach a target generalization gap."""
    if gap <= 0:
        raise ValueError("gap must be positive")
    t = max(int(n_trainable_gates), 1)
    numerator = t * np.log(t) if with_log and t > 1 else t
    return int(np.ceil(numerator / gap**2))

noise_survival

noise_survival(gate_fidelity: float, depth: int) -> float

f^d — independent per-gate success compounding over the circuit.

Source code in src/qmlkit/metrics.py

def noise_survival(gate_fidelity: float, depth: int) -> float:
    """``f^d`` — independent per-gate success compounding over the circuit."""
    return float(gate_fidelity**depth)

compare_ansatze

compare_ansatze(
    ansatze: Sequence[Ansatz],
    n_samples: int = 300,
    seed: int | None = 0,
) -> list[dict[str, object]]

The same report across candidates — the table a paper would want.

Source code in src/qmlkit/metrics.py

def compare_ansatze(
    ansatze: Sequence[Ansatz], n_samples: int = 300, seed: int | None = 0
) -> list[dict[str, object]]:
    """The same report across candidates — the table a paper would want."""
    return [AnsatzReport(a, n_samples=n_samples, seed=seed).results for a in ansatze]

qmlkit.fourier

fourier

What function does this model actually represent?

A variational model with data re-uploading is a truncated Fourier series in its inputs. The encoding fixes which frequencies are reachable; the ansatz only sets their coefficients. That is the central claim of the re-uploading literature, and this module turns it from a claim into a measurement:

coeffs = fourier_coefficients(f, degree=4)
spectrum(f)                 # which frequencies actually carry weight

Useful for two things. Checking that L uploads really did buy frequencies 0..L — and diagnosing a model that will not fit its target, because if the target's frequency is not in the spectrum, no amount of training will reach it.

fourier_coefficients

fourier_coefficients(
    f: ScalarFn,
    degree: int = 5,
    n_samples: int | None = None,
) -> NDArray[Any]

Complex Fourier coefficients c_-d .. c_d of a 2*pi-periodic function.

Sampled on a uniform grid and transformed exactly — no fitting, no optimiser. The grid must be at least 2*degree + 1 points to avoid aliasing, which is the default.

Source code in src/qmlkit/fourier.py

def fourier_coefficients(
    f: ScalarFn, degree: int = 5, n_samples: int | None = None
) -> npt.NDArray[Any]:
    """Complex Fourier coefficients ``c_-d .. c_d`` of a ``2*pi``-periodic function.

    Sampled on a uniform grid and transformed exactly — no fitting, no optimiser.
    The grid must be at least ``2*degree + 1`` points to avoid aliasing, which is
    the default.
    """
    if degree < 0:
        raise ValueError("degree must be non-negative")
    n = n_samples or (2 * degree + 1)
    if n < 2 * degree + 1:
        raise ValueError(
            f"{n} samples cannot resolve degree {degree}; need at least {2 * degree + 1}"
        )
    grid = np.linspace(0.0, 2 * np.pi, n, endpoint=False)
    values = np.array([float(f(x)) for x in grid], dtype=complex)
    full = np.fft.fft(values) / n
    return np.concatenate([full[-degree:], full[: degree + 1]]) if degree else full[:1]

spectrum

spectrum(
    f: ScalarFn,
    degree: int = 5,
    n_samples: int | None = None,
    tol: float = 1e-08,
) -> dict[int, float]

Frequency -> amplitude, keeping only what is actually present.

Amplitudes are |c_k| + |c_-k| for k > 0, so a real-valued model reports one number per frequency rather than a conjugate pair.

Source code in src/qmlkit/fourier.py

def spectrum(
    f: ScalarFn, degree: int = 5, n_samples: int | None = None, tol: float = 1e-8
) -> dict[int, float]:
    """Frequency -> amplitude, keeping only what is actually present.

    Amplitudes are ``|c_k| + |c_-k|`` for ``k > 0``, so a real-valued model reports
    one number per frequency rather than a conjugate pair.
    """
    coeffs = fourier_coefficients(f, degree, n_samples)
    out: dict[int, float] = {}
    zero = degree
    if abs(coeffs[zero]) > tol:
        out[0] = float(abs(coeffs[zero]))
    for k in range(1, degree + 1):
        amp = float(abs(coeffs[zero + k]) + abs(coeffs[zero - k]))
        if amp > tol:
            out[k] = amp
    return out

reconstruct

reconstruct(
    coeffs: NDArray[Any], x: float | NDArray[Any]
) -> NDArray[Any]

Evaluate the series these coefficients describe.

Source code in src/qmlkit/fourier.py

def reconstruct(coeffs: npt.NDArray[Any], x: float | npt.NDArray[Any]) -> npt.NDArray[Any]:
    """Evaluate the series these coefficients describe."""
    c = np.asarray(coeffs, dtype=complex)
    degree = (c.size - 1) // 2
    freqs = np.arange(-degree, degree + 1)
    xs = np.atleast_1d(np.asarray(x, dtype=float))
    return np.real(np.exp(1j * np.outer(xs, freqs)) @ c)

reachable_frequencies

reachable_frequencies(n_uploads: int) -> list[int]

L uploads of a Pauli-rotation encoding reach frequencies 0..L.

Source code in src/qmlkit/fourier.py

def reachable_frequencies(n_uploads: int) -> list[int]:
    """``L`` uploads of a Pauli-rotation encoding reach frequencies ``0..L``."""
    if n_uploads < 0:
        raise ValueError("n_uploads cannot be negative")
    return list(range(n_uploads + 1))

dominant_frequency

dominant_frequency(f: ScalarFn, degree: int = 8) -> int

The non-zero frequency carrying the most weight — 0 if the model is constant.

Source code in src/qmlkit/fourier.py

def dominant_frequency(f: ScalarFn, degree: int = 8) -> int:
    """The non-zero frequency carrying the most weight — 0 if the model is constant."""
    spec = {k: v for k, v in spectrum(f, degree).items() if k != 0}
    return max(spec, key=lambda k: spec[k]) if spec else 0

model_spectrum

model_spectrum(
    encoder: Any,
    theta: ArrayLike,
    obs: Observable | None = None,
    degree: int | None = None,
    backend: BackendLike = None,
) -> dict[int, float]

Spectrum of a one-feature re-uploading model, as a function of its input.

The direct check that an encoder buys the frequencies it claims to.

Source code in src/qmlkit/fourier.py

def model_spectrum(
    encoder: Any,
    theta: ArrayLike,
    obs: Observable | None = None,
    degree: int | None = None,
    backend: BackendLike = None,
) -> dict[int, float]:
    """Spectrum of a one-feature re-uploading model, as a function of its input.

    The direct check that an encoder buys the frequencies it claims to.
    """
    from qmlkit.core.execute import expval
    from qmlkit.core.observables import Z

    obs = Z(0) if obs is None else obs
    n_uploads = getattr(encoder, "n_uploads", 1)
    deg = degree if degree is not None else n_uploads + 2

    def f(x: float) -> float:
        return expval(encoder.build([x]), obs, theta=theta, backend=backend)

    return spectrum(f, deg)

qmlkit.info

info

Quantum information quantities — the qml.qinfo equivalent.

Reduced density matrices, purity, entropies and state fidelity. These are the building blocks the trainability metrics and the projected kernel are made of, and they are useful on their own for looking at what a circuit is actually doing.

density_matrix

density_matrix(
    state: CircuitSpec | NDArray[Any],
    backend: BackendLike = None,
) -> NDArray[Any]

|psi><psi| for a pure state.

Source code in src/qmlkit/info.py

def density_matrix(
    state: CircuitSpec | npt.NDArray[Any], backend: BackendLike = None
) -> npt.NDArray[Any]:
    """``|psi><psi|`` for a pure state."""
    psi = _as_state(state, backend)
    return np.outer(psi, psi.conj())

reduced_dm

reduced_dm(
    state: CircuitSpec | NDArray[Any],
    wires: Sequence[int],
    n_qubits: int | None = None,
    backend: BackendLike = None,
) -> NDArray[Any]

Trace out everything except wires, kept in the order given.

Qubit 0 is the most significant bit, matching the rest of the library, and the returned matrix carries its own bit order: qubit wires[0] is its most significant. reduced_dm(psi, [1, 0]) is therefore reduced_dm(psi, [0, 1]) with its two subsystems exchanged — SWAP @ rho @ SWAP — and not the same matrix. Everything basis-independent (trace, eigenvalues, purity, entropy) is blind to the difference, which is why the contract has to be stated here rather than left to be read off a test.

Source code in src/qmlkit/info.py

def reduced_dm(
    state: CircuitSpec | npt.NDArray[Any],
    wires: Sequence[int],
    n_qubits: int | None = None,
    backend: BackendLike = None,
) -> npt.NDArray[Any]:
    """Trace out everything except ``wires``, kept **in the order given**.

    Qubit 0 is the most significant bit, matching the rest of the library, and the
    returned matrix carries its own bit order: qubit ``wires[0]`` is its most
    significant. ``reduced_dm(psi, [1, 0])`` is therefore ``reduced_dm(psi, [0, 1])``
    with its two subsystems exchanged — ``SWAP @ rho @ SWAP`` — and not the same
    matrix. Everything basis-independent (trace, eigenvalues, purity, entropy) is
    blind to the difference, which is why the contract has to be stated here rather
    than left to be read off a test.
    """
    psi = _as_state(state, backend)
    n = n_qubits if n_qubits is not None else int(np.log2(psi.size))
    if 2**n != psi.size:
        raise ValueError(f"state of size {psi.size} is not {n} qubits")
    keep = [int(w) for w in wires]
    if any(not 0 <= w < n for w in keep):
        raise ValueError(f"wires {list(wires)} out of range for {n} qubits")
    twice = next((w for w in keep if keep.count(w) > 1), None)
    if twice is not None:
        raise ValueError(
            f"wires {list(wires)} names qubit {twice} twice: a subsystem holds each "
            "qubit once, so there is no reduced state this could mean"
        )

    tensor = psi.reshape((2,) * n)
    traced = [q for q in range(n) if q not in keep]
    # move kept axes to the front, then contract the rest against their conjugate
    perm = keep + traced
    tensor = np.transpose(tensor, perm)
    k = len(keep)
    tensor = tensor.reshape(2**k, -1)
    return tensor @ tensor.conj().T

purity

purity(
    state: CircuitSpec | NDArray[Any],
    wires: Sequence[int] | None = None,
    n_qubits: int | None = None,
    backend: BackendLike = None,
) -> float

Tr(rho^2) — 1 for a pure state, 1/d for the maximally mixed one.

With no wires the whole register is meant. That is 1 for a statevector by construction — but not when backend evolves a density matrix, where the whole point of asking is that the state may be mixed. Returning 1.0 there answered a question nobody asked, and was the one path in this module that neither refused nor computed.

Source code in src/qmlkit/info.py

def purity(
    state: CircuitSpec | npt.NDArray[Any],
    wires: Sequence[int] | None = None,
    n_qubits: int | None = None,
    backend: BackendLike = None,
) -> float:
    """``Tr(rho^2)`` — 1 for a pure state, ``1/d`` for the maximally mixed one.

    With no ``wires`` the whole register is meant. That is 1 for a statevector by
    construction — but not when ``backend`` evolves a density matrix, where the whole
    point of asking is that the state may be mixed. Returning 1.0 there answered a
    question nobody asked, and was the one path in this module that neither refused
    nor computed.
    """
    if wires is None:
        device = get_backend(backend)
        if not device.supports_statevector and hasattr(device, "purity"):
            if not isinstance(state, CircuitSpec):
                raise TypeError(
                    f"the {device.name!r} backend computes purity from a circuit, not "
                    "from an array; pass the CircuitSpec"
                )
            return float(device.purity(state))
        return 1.0  # a statevector is pure by construction
    rho = reduced_dm(state, wires, n_qubits, backend)
    return float(np.real(np.trace(rho @ rho)))

vn_entropy

vn_entropy(
    state: CircuitSpec | NDArray[Any],
    wires: Sequence[int],
    n_qubits: int | None = None,
    base: float | None = None,
    backend: BackendLike = None,
) -> float

Von Neumann entropy of a subsystem — how entangled it is with the rest.

Source code in src/qmlkit/info.py

def vn_entropy(
    state: CircuitSpec | npt.NDArray[Any],
    wires: Sequence[int],
    n_qubits: int | None = None,
    base: float | None = None,
    backend: BackendLike = None,
) -> float:
    """Von Neumann entropy of a subsystem — how entangled it is with the rest."""
    rho = reduced_dm(state, wires, n_qubits, backend)
    eig = np.linalg.eigvalsh(rho)
    eig = eig[eig > 1e-12]
    entropy = float(-np.sum(eig * np.log(eig)))
    return entropy / np.log(base) if base is not None else entropy

mutual_info

mutual_info(
    state: CircuitSpec | NDArray[Any],
    wires_a: Sequence[int],
    wires_b: Sequence[int],
    n_qubits: int | None = None,
    backend: BackendLike = None,
) -> float

S(A) + S(B) - S(AB) — total correlation between two subsystems.

Source code in src/qmlkit/info.py

def mutual_info(
    state: CircuitSpec | npt.NDArray[Any],
    wires_a: Sequence[int],
    wires_b: Sequence[int],
    n_qubits: int | None = None,
    backend: BackendLike = None,
) -> float:
    """``S(A) + S(B) - S(AB)`` — total correlation between two subsystems."""
    a, b = sorted(set(wires_a)), sorted(set(wires_b))
    if set(a) & set(b):
        raise ValueError(f"subsystems overlap on {sorted(set(a) & set(b))}")
    sa = vn_entropy(state, a, n_qubits, backend=backend)
    sb = vn_entropy(state, b, n_qubits, backend=backend)
    sab = vn_entropy(state, a + b, n_qubits, backend=backend)
    return sa + sb - sab

state_fidelity

state_fidelity(
    state_a: CircuitSpec | NDArray[Any],
    state_b: CircuitSpec | NDArray[Any],
    backend: BackendLike = None,
) -> float

|<a|b>|^2 — the quantity a fidelity kernel estimates.

Source code in src/qmlkit/info.py

def state_fidelity(
    state_a: CircuitSpec | npt.NDArray[Any],
    state_b: CircuitSpec | npt.NDArray[Any],
    backend: BackendLike = None,
) -> float:
    """``|<a|b>|^2`` — the quantity a fidelity kernel estimates."""
    a = _as_state(state_a, backend)
    b = _as_state(state_b, backend)
    if a.size != b.size:
        raise ValueError(f"states have different widths: {a.size} vs {b.size}")
    return float(abs(np.vdot(a, b)) ** 2)

concurrence

concurrence(
    state: CircuitSpec | NDArray[Any],
    backend: BackendLike = None,
) -> float

Two-qubit concurrence — 0 for a product state, 1 for a Bell state.

Source code in src/qmlkit/info.py

def concurrence(state: CircuitSpec | npt.NDArray[Any], backend: BackendLike = None) -> float:
    """Two-qubit concurrence — 0 for a product state, 1 for a Bell state."""
    psi = _as_state(state, backend)
    if psi.size != 4:
        raise ValueError("concurrence is defined here for two qubits only")
    yy = np.array([[0, 0, 0, -1], [0, 0, 1, 0], [0, 1, 0, 0], [-1, 0, 0, 0]], dtype=complex)
    return float(abs(psi @ yy @ psi))

bloch_vector

bloch_vector(
    state: CircuitSpec | NDArray[Any],
    wire: int = 0,
    n_qubits: int | None = None,
    backend: BackendLike = None,
) -> NDArray[Any]

(<X>, <Y>, <Z>) for one qubit — its point on (or in) the Bloch sphere.

Source code in src/qmlkit/info.py

def bloch_vector(
    state: CircuitSpec | npt.NDArray[Any],
    wire: int = 0,
    n_qubits: int | None = None,
    backend: BackendLike = None,
) -> npt.NDArray[Any]:
    """``(<X>, <Y>, <Z>)`` for one qubit — its point on (or in) the Bloch sphere."""
    rho = reduced_dm(state, [wire], n_qubits, backend)
    x = 2 * np.real(rho[0, 1])
    y = 2 * np.imag(rho[1, 0])
    z = np.real(rho[0, 0] - rho[1, 1])
    return np.array([x, y, z], dtype=float)

qmlkit.optim

optim

Optimisers that only make sense for quantum circuits.

Adam and SGD come from torch. These three do not exist there, because they exploit structure a general optimiser cannot see:

• Rotosolve — a circuit expectation is a sinusoid in any single Pauli-rotation angle. Three evaluations pin that sinusoid down exactly, so you can jump straight to its minimum instead of stepping toward it. No learning rate, no tuning.
• Quantum natural gradient — parameter space is curved. Following the Fubini–Study geometry rather than the Euclidean one usually converges in far fewer steps.
• SPSA — lives in :mod:qmlkit.gradients.spsa; two evaluations per step at any parameter count.

AdamState dataclass

AdamState(m: NDArray[Any], v: NDArray[Any], t: int = 0)

The two running averages Adam carries between steps, and the step count.

Exposed because a circuit-level optimiser is often driven one step at a time from a loop the caller owns, and losing this between steps silently turns Adam back into plain gradient descent with a decaying learning rate.

supports_rotosolve

supports_rotosolve(spec: CircuitSpec) -> bool

Whether Rotosolve's closed form is actually valid for this circuit.

Rotosolve works because a circuit expectation is a single sinusoid :math:A\sin(\theta + B) + C in any one Pauli-rotation angle — three samples then determine it exactly. That holds when a parameter drives one rotation, or several that compose into one (same qubit, same generator).

It does not hold when one parameter is shared across gates that do not compose — QAOA's cost angle drives one rz per graph edge, so E(gamma) carries one frequency per edge. Rotosolve then solves for the wrong minimum, converges immediately, and reports a number that looks like a result. Measured on a 5-edge MaxCut: frequencies 1 through 4 are all present, and Rotosolve sticks at the uniform-state energy no matter how many sweeps it is given.

Source code in src/qmlkit/optim.py

def supports_rotosolve(spec: CircuitSpec) -> bool:
    r"""Whether Rotosolve's closed form is actually valid for this circuit.

    Rotosolve works because a circuit expectation is a *single* sinusoid
    :math:`A\sin(\theta + B) + C` in any one Pauli-rotation angle — three samples then
    determine it exactly. That holds when a parameter drives one rotation, or several
    that compose into one (same qubit, same generator).

    It does **not** hold when one parameter is shared across gates that do not compose
    — QAOA's cost angle drives one ``rz`` per graph edge, so ``E(gamma)`` carries one
    frequency per edge. Rotosolve then solves for the wrong minimum, converges
    immediately, and reports a number that looks like a result. Measured on a 5-edge
    MaxCut: frequencies 1 through 4 are all present, and Rotosolve sticks at the
    uniform-state energy no matter how many sweeps it is given.
    """
    from collections import defaultdict

    by_param: dict[int, list[tuple[str, tuple[int, ...]]]] = defaultdict(list)
    for slot in spec.slots():
        op = spec.ops[slot.op_index]
        by_param[slot.ref.index].append((slot.gate, op.qubits))
    # several occurrences are fine only if they are the same rotation on the same wire
    return all(len(set(sites)) <= 1 for sites in by_param.values())

rotosolve_step

rotosolve_step(
    f: LossFn,
    theta: ArrayLike,
    indices: Sequence[int] | None = None,
) -> NDArray[Any]

One sweep: set every coordinate to its exact optimum, in turn.

Three evaluations per parameter, and each one lands on that coordinate's minimum rather than moving toward it.

Source code in src/qmlkit/optim.py

def rotosolve_step(
    f: LossFn, theta: ArrayLike, indices: Sequence[int] | None = None
) -> npt.NDArray[Any]:
    """One sweep: set every coordinate to its exact optimum, in turn.

    Three evaluations per parameter, and each one lands on that coordinate's minimum
    rather than moving toward it.
    """
    arr = np.asarray(theta, dtype=float).ravel().copy()
    for k in indices if indices is not None else range(arr.size):
        arr[k] = _optimal_angle(f, arr, k)
    return arr

minimize_rotosolve

minimize_rotosolve(
    f: LossFn,
    theta0: Sequence[float],
    n_sweeps: int = 20,
    tol: float = 1e-09,
    callback: Callable[[int, NDArray[Any], float], None]
    | None = None,
) -> tuple[NDArray[Any], list[float]]

Minimise by repeated Rotosolve sweeps. No learning rate to choose.

Precondition. f must be a single sinusoid in each angle — true for a plain expectation value <O> of a circuit where each parameter drives one Pauli rotation. It is not true when an angle is shared across gates that do not compose (QAOA's cost angle drives one rz per edge), nor when the loss is non-linear in the state (purity is Tr(rho^2), so it carries double frequencies). In those cases Rotosolve converges immediately on the wrong point and reports it as a result. :func:supports_rotosolve checks the first case; the second is a property of your loss, not of the circuit.

Source code in src/qmlkit/optim.py

def minimize_rotosolve(
    f: LossFn,
    theta0: Sequence[float],
    n_sweeps: int = 20,
    tol: float = 1e-9,
    callback: Callable[[int, npt.NDArray[Any], float], None] | None = None,
) -> tuple[npt.NDArray[Any], list[float]]:
    """Minimise by repeated Rotosolve sweeps. No learning rate to choose.

    **Precondition.** ``f`` must be a single sinusoid in each angle — true for a plain
    expectation value ``<O>`` of a circuit where each parameter drives one Pauli
    rotation. It is *not* true when an angle is shared across gates that do not
    compose (QAOA's cost angle drives one ``rz`` per edge), nor when the loss is
    non-linear in the state (purity is ``Tr(rho^2)``, so it carries double
    frequencies). In those cases Rotosolve converges immediately on the wrong point
    and reports it as a result. :func:`supports_rotosolve` checks the first case;
    the second is a property of your loss, not of the circuit.
    """
    theta: npt.NDArray[Any] = np.asarray(theta0, dtype=float).ravel().copy()
    history = [float(f(theta))]
    for sweep in range(n_sweeps):
        theta = rotosolve_step(f, theta)
        value = float(f(theta))
        history.append(value)
        if callback is not None:
            callback(sweep, theta, value)
        if abs(history[-2] - value) < tol:
            break
    return theta, history

adam_step

adam_step(
    theta: ArrayLike,
    gradient: ArrayLike,
    state: AdamState,
    lr: float = 0.05,
    beta1: float = 0.9,
    beta2: float = 0.999,
    eps: float = 1e-08,
) -> tuple[NDArray[Any], AdamState]

One Adam update, given a gradient you already have.

Returns the new parameters and the new state; state is not mutated, so a caller can keep a trajectory without copying by hand.

Adam earns its place on variational circuits for a specific reason: parameter gradients in a deep ansatz differ in scale by orders of magnitude — a rotation near the readout moves the expectation far more than one behind a wall of entanglers — and a single learning rate either crawls on the small ones or diverges on the large. Dividing by the running gradient magnitude makes the step size per-parameter, which is exactly that problem.

Source code in src/qmlkit/optim.py

def adam_step(
    theta: ArrayLike,
    gradient: ArrayLike,
    state: AdamState,
    lr: float = 0.05,
    beta1: float = 0.9,
    beta2: float = 0.999,
    eps: float = 1e-8,
) -> tuple[npt.NDArray[Any], AdamState]:
    """One Adam update, given a gradient you already have.

    Returns the new parameters and the new state; ``state`` is not mutated, so a
    caller can keep a trajectory without copying by hand.

    Adam earns its place on variational circuits for a specific reason: parameter
    gradients in a deep ansatz differ in scale by orders of magnitude — a rotation
    near the readout moves the expectation far more than one behind a wall of
    entanglers — and a single learning rate either crawls on the small ones or
    diverges on the large. Dividing by the running gradient magnitude makes the step
    size per-parameter, which is exactly that problem.
    """
    values = np.asarray(theta, dtype=float).ravel()
    grad = np.asarray(gradient, dtype=float).ravel()
    if grad.shape != values.shape:
        raise ValueError(f"gradient has {grad.size} entries but there are {values.size} parameters")

    t = state.t + 1
    m = beta1 * state.m + (1.0 - beta1) * grad
    v = beta2 * state.v + (1.0 - beta2) * grad**2
    # bias correction: m and v start at zero, so without this the first steps are
    # damped by roughly (1 - beta) and Adam looks like it is not moving
    m_hat = m / (1.0 - beta1**t)
    v_hat = v / (1.0 - beta2**t)
    stepped: npt.NDArray[Any] = values - lr * m_hat / (np.sqrt(v_hat) + eps)
    return stepped, AdamState(m=m, v=v, t=t)

minimize_adam

minimize_adam(
    f: LossFn,
    theta0: Sequence[float],
    grad: Callable[[NDArray[Any]], NDArray[Any]],
    n_steps: int = 100,
    lr: float = 0.05,
    beta1: float = 0.9,
    beta2: float = 0.999,
    eps: float = 1e-08,
    tol: float = 0.0,
    callback: Callable[[int, NDArray[Any], float], None]
    | None = None,
) -> tuple[NDArray[Any], list[float]]

Minimise f by Adam, using the gradient grad supplies.

grad is explicit rather than inferred because the right way to differentiate a circuit depends on the circuit and the backend: pass lambda t: qk.grad(spec, t, obs) for an exact gradient, or a shot-based one when the point is to see what a device would do. :func:~qmlkit.grad chooses the method; this chooses the step.

tol stops early when the loss improves by less than that between steps; the default of 0.0 runs the full budget, since a variational loss plateaus and then escapes often enough that stopping on the first flat step is usually wrong.

Source code in src/qmlkit/optim.py

def minimize_adam(
    f: LossFn,
    theta0: Sequence[float],
    grad: Callable[[npt.NDArray[Any]], npt.NDArray[Any]],
    n_steps: int = 100,
    lr: float = 0.05,
    beta1: float = 0.9,
    beta2: float = 0.999,
    eps: float = 1e-8,
    tol: float = 0.0,
    callback: Callable[[int, npt.NDArray[Any], float], None] | None = None,
) -> tuple[npt.NDArray[Any], list[float]]:
    """Minimise ``f`` by Adam, using the gradient ``grad`` supplies.

    ``grad`` is explicit rather than inferred because the right way to differentiate a
    circuit depends on the circuit and the backend: pass
    ``lambda t: qk.grad(spec, t, obs)`` for an exact gradient, or a shot-based one
    when the point is to see what a device would do. :func:`~qmlkit.grad` chooses the
    method; this chooses the step.

    ``tol`` stops early when the loss improves by less than that between steps;
    the default of ``0.0`` runs the full budget, since a variational loss plateaus and
    then escapes often enough that stopping on the first flat step is usually wrong.
    """
    theta: npt.NDArray[Any] = np.asarray(theta0, dtype=float).ravel().copy()
    state = AdamState.for_parameters(theta.size)
    history = [float(f(theta))]
    for step in range(n_steps):
        theta, state = adam_step(theta, grad(theta), state, lr, beta1, beta2, eps)
        value = float(f(theta))
        history.append(value)
        if callback is not None:
            callback(step, theta, value)
        if tol > 0.0 and abs(history[-2] - value) < tol:
            break
    return theta, history

metric_tensor

metric_tensor(
    spec: CircuitSpec,
    theta: ArrayLike,
    approx: str | None = "block-diag",
    backend: BackendLike = None,
    eps: float = 0.0001,
) -> NDArray[Any]

Fubini–Study metric — the curvature of parameter space.

approx="diag" keeps only the diagonal (cheapest). "block-diag" and None compute the full tensor from state overlaps; on a simulator that is affordable and exact, so they currently coincide. Note that PennyLane's approx="block-diag" means something narrower — it blocks by circuit layer and zeroes every cross-layer entry — so the same keyword does not port between the two libraries. qmlkit follows the true geometry; PennyLane follows an approximation to it.

.. math:: g_{ij} = \mathrm{Re}\langle \partial_i\psi | \partial_j\psi \rangle - \langle \partial_i\psi|\psi\rangle\langle\psi|\partial_j\psi\rangle

The derivative states are exact whenever every parameterised gate declares a closed-form derivative, which every built-in gate does. eps is only consulted on the fallback path, for a custom gate registered without a dmatrix.

Source code in src/qmlkit/optim.py

def metric_tensor(
    spec: CircuitSpec,
    theta: ArrayLike,
    approx: str | None = "block-diag",
    backend: BackendLike = None,
    eps: float = 1e-4,
) -> npt.NDArray[Any]:
    r"""Fubini–Study metric — the curvature of parameter space.

    ``approx="diag"`` keeps only the diagonal (cheapest). ``"block-diag"`` and
    ``None`` compute the full tensor from state overlaps; on a simulator that is
    affordable and exact, so they currently coincide. Note that PennyLane's
    ``approx="block-diag"`` means something narrower — it blocks by circuit layer and
    zeroes every cross-layer entry — so the same keyword does not port between the
    two libraries. qmlkit follows the true geometry; PennyLane follows an
    approximation to it.

    .. math::  g_{ij} = \mathrm{Re}\langle \partial_i\psi | \partial_j\psi \rangle
               - \langle \partial_i\psi|\psi\rangle\langle\psi|\partial_j\psi\rangle

    The derivative states are exact whenever every parameterised gate declares a
    closed-form derivative, which every built-in gate does. ``eps`` is only consulted
    on the fallback path, for a custom gate registered without a ``dmatrix``.
    """
    from qmlkit.core.execute import statevector
    from qmlkit.gradients.adjoint import supports_adjoint

    require_statevector(backend, "the Fubini-Study metric tensor")
    arr = np.asarray(theta, dtype=float).ravel()
    p = arr.size

    if supports_adjoint(spec, backend):
        psi, derivatives = _exact_derivative_states(spec, arr, backend)
    else:
        psi = statevector(spec.bind(arr), backend=backend)
        derivatives = np.empty((p, psi.size), dtype=complex)
        for k in range(p):
            plus, minus = arr.copy(), arr.copy()
            plus[k] += eps
            minus[k] -= eps
            derivatives[k] = (
                statevector(spec.bind(plus), backend=backend)
                - statevector(spec.bind(minus), backend=backend)
            ) / (2 * eps)

    overlaps = derivatives @ psi.conj()
    g = np.real(derivatives.conj() @ derivatives.T) - np.real(np.outer(overlaps.conj(), overlaps))
    if approx == "diag":
        return np.diag(np.diag(g))
    if approx in ("block-diag", None):
        return g
    raise unknown(
        "approx", approx, ("diag", "block-diag"), hint="Pass None for the full metric tensor."
    )

quantum_fisher_information

quantum_fisher_information(
    spec: CircuitSpec,
    theta: Sequence[float],
    backend: BackendLike = None,
) -> NDArray[Any]

QFIM — exactly 4 x the Fubini–Study metric.

Distinct from the classical Fisher information in :mod:qmlkit.metrics, which describes the output distribution and is what effective dimension uses.

Source code in src/qmlkit/optim.py

def quantum_fisher_information(
    spec: CircuitSpec, theta: Sequence[float], backend: BackendLike = None
) -> npt.NDArray[Any]:
    """QFIM — exactly ``4 x`` the Fubini–Study metric.

    Distinct from the *classical* Fisher information in :mod:`qmlkit.metrics`, which
    describes the output distribution and is what effective dimension uses.
    """
    return 4.0 * metric_tensor(spec, theta, None, backend)

qng_step

qng_step(
    spec: CircuitSpec,
    theta: ArrayLike,
    obs: Observable | None = None,
    lr: float = 0.1,
    approx: str = "block-diag",
    regularization: float = 1e-06,
    backend: BackendLike = None,
) -> NDArray[Any]

One natural-gradient step: theta <- theta - lr * g^+ grad.

The pseudo-inverse of the metric rescales each direction by how much the state actually moves, rather than how much the parameter does.

Source code in src/qmlkit/optim.py

def qng_step(
    spec: CircuitSpec,
    theta: ArrayLike,
    obs: Observable | None = None,
    lr: float = 0.1,
    approx: str = "block-diag",
    regularization: float = 1e-6,
    backend: BackendLike = None,
) -> npt.NDArray[Any]:
    r"""One natural-gradient step: ``theta <- theta - lr * g^+ grad``.

    The pseudo-inverse of the metric rescales each direction by how much the *state*
    actually moves, rather than how much the parameter does.
    """
    from qmlkit.gradients.dispatch import grad

    obs = Z(0) if obs is None else obs
    arr = np.asarray(theta, dtype=float).ravel()
    g = metric_tensor(spec, arr, approx, backend)
    gradient = grad(spec, arr, obs, backend=backend)
    natural = np.linalg.pinv(g + regularization * np.eye(g.shape[0])) @ gradient
    return arr - lr * natural

minimize_qng

minimize_qng(
    spec: CircuitSpec,
    theta0: Sequence[float],
    obs: Observable | None = None,
    n_steps: int = 50,
    lr: float = 0.1,
    approx: str = "block-diag",
    backend: BackendLike = None,
    callback: Callable[[int, NDArray[Any], float], None]
    | None = None,
) -> tuple[NDArray[Any], list[float]]

Minimise <obs> by quantum natural gradient descent.

Source code in src/qmlkit/optim.py

def minimize_qng(
    spec: CircuitSpec,
    theta0: Sequence[float],
    obs: Observable | None = None,
    n_steps: int = 50,
    lr: float = 0.1,
    approx: str = "block-diag",
    backend: BackendLike = None,
    callback: Callable[[int, npt.NDArray[Any], float], None] | None = None,
) -> tuple[npt.NDArray[Any], list[float]]:
    """Minimise ``<obs>`` by quantum natural gradient descent."""
    obs = Z(0) if obs is None else obs
    theta: npt.NDArray[Any] = np.asarray(theta0, dtype=float).ravel().copy()
    history: list[float] = []
    for step in range(n_steps):
        value = expval(spec, obs, theta=theta, backend=backend)
        history.append(value)
        if callback is not None:
            callback(step, theta, value)
        theta = qng_step(spec, theta, obs, lr, approx, backend=backend)
    history.append(expval(spec, obs, theta=theta, backend=backend))
    return theta, history

shots_for_precision

shots_for_precision(eps: float, z: float = 0.0) -> int

Shots needed to reach standard error eps — the 1/eps**2 price.

Source code in src/qmlkit/utils/shots.py

def shots_for_precision(eps: float, z: float = 0.0) -> int:
    """Shots needed to reach standard error ``eps`` — the ``1/eps**2`` price."""
    if eps <= 0:
        raise ValueError("eps must be positive")
    return int(np.ceil(variance(z) / eps**2))

qmlkit.datasets

datasets

Datasets for benchmarking quantum models.

Small, self-contained, no downloads, no sklearn. The important one is :func:ad_hoc_data, which is constructed to be separable by a specific quantum feature map and not by a classical kernel — so it distinguishes a working implementation from one that only appears to work.

ad_hoc_data

ad_hoc_data(
    n_samples: int = 40,
    n_features: int = 2,
    gap: float = 0.3,
    seed: int | None = None,
    scale: float = 2 * pi,
) -> tuple[NDArray[Any], NDArray[Any]]

The Havlíček-style separable-by-construction dataset.

Labels come from the sign of a hidden observable measured on a ZZ-feature-mapped state, so the ZZ kernel separates it by construction while classical kernels struggle. gap discards points near the boundary, which makes the separation clean enough to be a real check.

Source code in src/qmlkit/datasets.py

def ad_hoc_data(
    n_samples: int = 40,
    n_features: int = 2,
    gap: float = 0.3,
    seed: int | None = None,
    scale: float = 2 * np.pi,
) -> tuple[npt.NDArray[Any], npt.NDArray[Any]]:
    """The Havlíček-style separable-by-construction dataset.

    Labels come from the sign of a hidden observable measured on a ZZ-feature-mapped
    state, so the ZZ kernel separates it by construction while classical kernels
    struggle. ``gap`` discards points near the boundary, which makes the separation
    clean enough to be a real check.
    """
    from qmlkit.core.execute import expectation
    from qmlkit.core.observables import PauliString
    from qmlkit.encoding.feature_maps import ZZFeatureMap

    rng = np.random.default_rng(seed)
    fmap = ZZFeatureMap(n_features, reps=2)
    witness = PauliString(tuple((q, "Z") for q in range(n_features)))

    xs: list[npt.NDArray[Any]] = []
    ys: list[int] = []
    attempts = 0
    while len(xs) < n_samples and attempts < 200 * n_samples:
        attempts += 1
        x = rng.uniform(0, scale, n_features)
        value = float(expectation(fmap.build(x), witness))
        if abs(value) < gap:
            continue  # too close to the boundary to label cleanly
        xs.append(x)
        ys.append(1 if value > 0 else 0)
    if len(xs) < n_samples:  # pragma: no cover - only with an extreme gap
        raise ValueError(f"could not find {n_samples} samples outside a gap of {gap}")
    return np.array(xs), np.array(ys)

bars_and_stripes

bars_and_stripes(
    size: int = 2, seed: int | None = None
) -> NDArray[Any]

Every bars-and-stripes pattern on a size x size grid, flattened.

The standard target distribution for a quantum circuit Born machine: a sparse, highly structured subset of all bitstrings.

Source code in src/qmlkit/datasets.py

def bars_and_stripes(size: int = 2, seed: int | None = None) -> npt.NDArray[Any]:
    """Every bars-and-stripes pattern on a ``size x size`` grid, flattened.

    The standard target distribution for a quantum circuit Born machine: a sparse,
    highly structured subset of all bitstrings.
    """
    if size < 1:
        raise ValueError("size must be at least 1")
    patterns = set()
    for bits in itertools.product([0, 1], repeat=size):
        grid = np.tile(np.array(bits)[:, None], (1, size))  # bars
        patterns.add(tuple(grid.ravel()))
        patterns.add(tuple(grid.T.ravel()))  # stripes
    out = np.array(sorted(patterns), dtype=int)
    if seed is not None:
        np.random.default_rng(seed).shuffle(out)
    return out

make_moons

make_moons(
    n_samples: int = 100,
    noise: float = 0.1,
    seed: int | None = None,
    to_angles: bool = True,
) -> tuple[NDArray[Any], NDArray[Any]]

Two interleaving half-circles — not linearly separable.

Source code in src/qmlkit/datasets.py

def make_moons(
    n_samples: int = 100, noise: float = 0.1, seed: int | None = None, to_angles: bool = True
) -> tuple[npt.NDArray[Any], npt.NDArray[Any]]:
    """Two interleaving half-circles — not linearly separable."""
    rng = np.random.default_rng(seed)
    n_out = n_samples // 2
    n_in = n_samples - n_out
    t_out = np.linspace(0, np.pi, n_out)
    t_in = np.linspace(0, np.pi, n_in)
    x = np.vstack(
        [
            np.column_stack([np.cos(t_out), np.sin(t_out)]),
            np.column_stack([1 - np.cos(t_in), 1 - np.sin(t_in) - 0.5]),
        ]
    )
    x += rng.normal(0, noise, x.shape)
    y = np.array([0] * n_out + [1] * n_in)
    return (_to_angles(x) if to_angles else x), y

make_circles

make_circles(
    n_samples: int = 100,
    noise: float = 0.08,
    factor: float = 0.5,
    seed: int | None = None,
    to_angles: bool = True,
) -> tuple[NDArray[Any], NDArray[Any]]

One circle inside another — needs a nonlinear boundary.

Source code in src/qmlkit/datasets.py

def make_circles(
    n_samples: int = 100,
    noise: float = 0.08,
    factor: float = 0.5,
    seed: int | None = None,
    to_angles: bool = True,
) -> tuple[npt.NDArray[Any], npt.NDArray[Any]]:
    """One circle inside another — needs a nonlinear boundary."""
    rng = np.random.default_rng(seed)
    n_out = n_samples // 2
    n_in = n_samples - n_out
    t_out = np.linspace(0, 2 * np.pi, n_out, endpoint=False)
    t_in = np.linspace(0, 2 * np.pi, n_in, endpoint=False)
    x = np.vstack(
        [
            np.column_stack([np.cos(t_out), np.sin(t_out)]),
            factor * np.column_stack([np.cos(t_in), np.sin(t_in)]),
        ]
    )
    x += rng.normal(0, noise, x.shape)
    y = np.array([0] * n_out + [1] * n_in)
    return (_to_angles(x) if to_angles else x), y

make_blobs

make_blobs(
    n_samples: int = 100,
    centers: int = 2,
    spread: float = 0.4,
    n_features: int = 2,
    seed: int | None = None,
    to_angles: bool = True,
) -> tuple[NDArray[Any], NDArray[Any]]

Gaussian clusters — the easy baseline every model should pass.

Source code in src/qmlkit/datasets.py

def make_blobs(
    n_samples: int = 100,
    centers: int = 2,
    spread: float = 0.4,
    n_features: int = 2,
    seed: int | None = None,
    to_angles: bool = True,
) -> tuple[npt.NDArray[Any], npt.NDArray[Any]]:
    """Gaussian clusters — the easy baseline every model should pass."""
    rng = np.random.default_rng(seed)
    middles = rng.uniform(-2, 2, (centers, n_features))
    per = n_samples // centers
    xs, ys = [], []
    for c in range(centers):
        count = per if c < centers - 1 else n_samples - per * (centers - 1)
        xs.append(rng.normal(middles[c], spread, (count, n_features)))
        ys.extend([c] * count)
    x = np.vstack(xs)
    return (_to_angles(x) if to_angles else x), np.array(ys)

make_parity

make_parity(
    n_samples: int = 100,
    n_features: int = 4,
    seed: int | None = None,
) -> tuple[NDArray[Any], NDArray[Any]]

Label is the parity of the bits — the classic hard case for shallow models.

Source code in src/qmlkit/datasets.py

def make_parity(
    n_samples: int = 100, n_features: int = 4, seed: int | None = None
) -> tuple[npt.NDArray[Any], npt.NDArray[Any]]:
    """Label is the parity of the bits — the classic hard case for shallow models."""
    rng = np.random.default_rng(seed)
    bits = rng.integers(0, 2, (n_samples, n_features))
    y = bits.sum(axis=1) % 2
    return bits * np.pi, y  # 0 or pi, already an angle

train_test_split

train_test_split(
    X: NDArray[Any],
    y: NDArray[Any],
    test_size: float = 0.3,
    seed: int | None = None,
) -> tuple[
    NDArray[Any], NDArray[Any], NDArray[Any], NDArray[Any]
]

Shuffle and split. Here so a quickstart needs no extra dependency.

Source code in src/qmlkit/datasets.py

def train_test_split(
    X: npt.NDArray[Any], y: npt.NDArray[Any], test_size: float = 0.3, seed: int | None = None
) -> tuple[npt.NDArray[Any], npt.NDArray[Any], npt.NDArray[Any], npt.NDArray[Any]]:
    """Shuffle and split. Here so a quickstart needs no extra dependency."""
    rows = np.atleast_2d(np.asarray(X))
    labels = np.asarray(y).ravel()
    if not 0 < test_size < 1:
        raise ValueError("test_size must be between 0 and 1")
    rng = np.random.default_rng(seed)
    order = rng.permutation(len(rows))
    cut = int(round(len(rows) * (1 - test_size)))
    tr, te = order[:cut], order[cut:]
    return rows[tr], rows[te], labels[tr], labels[te]

qmlkit.draw

draw

Look at the circuit — the qml.draw / qml.specs equivalent.

Plain text, no matplotlib, no optional dependency. Reads the IR, so it works for any circuit the library can build and shows exactly what a backend will run.

draw

draw(
    spec: CircuitSpec,
    max_width: int = 160,
    ascii: bool | None = None,
) -> str

A text diagram of the circuit.

print(qk.draw(qk.hardware_efficient(3, 1).build()))

The diagram uses box-drawing glyphs. ascii=None (the default) checks whether sys.stdout can encode them and degrades to -, | and t when it cannot: a Windows console defaults to cp1252, which encodes none of them, and printing the result would otherwise raise UnicodeEncodeError from inside the caller, with a traceback naming the codec rather than this function. Pass ascii=True or ascii=False to decide for yourself. Column widths are identical either way, so the fallback lines up exactly like the Unicode form.

Source code in src/qmlkit/draw.py

def draw(spec: CircuitSpec, max_width: int = 160, ascii: bool | None = None) -> str:
    """A text diagram of the circuit.

    print(qk.draw(qk.hardware_efficient(3, 1).build()))

    The diagram uses box-drawing glyphs. ``ascii=None`` (the default) checks whether
    ``sys.stdout`` can encode them and degrades to ``-``, ``|`` and ``t`` when it
    cannot: a Windows console defaults to cp1252, which encodes none of them, and
    printing the result would otherwise raise ``UnicodeEncodeError`` from inside the
    caller, with a traceback naming the codec rather than this function. Pass
    ``ascii=True`` or ``ascii=False`` to decide for yourself. Column widths are
    identical either way, so the fallback lines up exactly like the Unicode form.
    """
    n = spec.n_qubits
    # pack operations into columns so nothing overlaps on a wire
    columns: list[dict[int, str]] = []
    frontier = [0] * n
    spans: list[tuple[int, int, int]] = []  # (column, top, bottom) for two-qubit links

    for op in spec.ops:
        col = max(frontier[q] for q in op.qubits)
        while len(columns) <= col:
            columns.append({})
        g = get_gate(op.gate)
        if g.n_qubits == 1:
            columns[col][op.qubits[0]] = _cell(op.gate, op.params)
        else:
            control, target = op.qubits[0], op.qubits[1]
            columns[col][control] = "@" if op.gate != "swap" else "x"
            columns[col][target] = _cell(_LABEL.get(op.gate, op.gate), op.params)
            spans.append((col, min(op.qubits), max(op.qubits)))
        for q in range(min(op.qubits), max(op.qubits) + 1):
            frontier[q] = col + 1

    widths = [max((len(v) for v in col.values()), default=1) for col in columns] or [1]
    label_w = len(f"q{n - 1}: ")

    wires = [f"q{q}: ".ljust(label_w) for q in range(n)]
    links = {(c, q) for c, top, bottom in spans for q in range(top, bottom + 1)}

    # `cells`, not `col`: `col` is a column *index* in the packing loop above, and
    # reusing the name for the column's contents is what made this look ill-typed
    for c, cells in enumerate(columns):
        w = widths[c]
        for q in range(n):
            if q in cells:
                wires[q] += "─" + cells[q].center(w, "─") + "─"
            elif (c, q) in links:
                wires[q] += "─" + "│".center(w, "─") + "─"
            else:
                wires[q] += "─" * (w + 2)

    out = [w + "─" for w in wires]
    if max(len(line) for line in out) > max_width:
        out = [line[: max_width - 3] + "..." for line in out]
    diagram = "\n".join(out)
    plain = not _stream_handles_unicode() if ascii is None else ascii
    return _downgrade(diagram) if plain else diagram

specs

specs(spec: CircuitSpec) -> dict[str, object]

Everything worth knowing about a circuit's cost, in one dict.

Source code in src/qmlkit/draw.py

def specs(spec: CircuitSpec) -> dict[str, object]:
    """Everything worth knowing about a circuit's cost, in one dict."""
    from qmlkit.gradients.parameter_shift import grad_circuit_cost

    out = dict(spec.resources())
    out["grad_circuits_parameter_shift"] = grad_circuit_cost(spec)
    out["grad_passes_adjoint"] = 1
    occurrences = {i: len(spec.occurrences_of(i)) for i in range(spec.n_params)}
    out["n_occurrences"] = occurrences
    tied = {i: c for i, c in occurrences.items() if c > 1}
    out["weight_tied_parameters"] = len(tied)
    return out

draw_ansatz

draw_ansatz(ansatz: object, max_width: int = 160) -> str

Convenience: draw an :class:~qmlkit.ansatz.library.Ansatz unbound.

Source code in src/qmlkit/draw.py

def draw_ansatz(ansatz: object, max_width: int = 160) -> str:
    """Convenience: draw an :class:`~qmlkit.ansatz.library.Ansatz` unbound."""
    return draw(ansatz.build(), max_width)  # type: ignore[attr-defined]

probabilities_bar

probabilities_bar(
    probs: NDArray[Any],
    n_qubits: int,
    top: int = 8,
    width: int = 30,
    ascii: bool | None = None,
) -> str

A text histogram of outcome probabilities — the most likely bitstrings.

Source code in src/qmlkit/draw.py

def probabilities_bar(
    probs: npt.NDArray[Any],
    n_qubits: int,
    top: int = 8,
    width: int = 30,
    ascii: bool | None = None,
) -> str:
    """A text histogram of outcome probabilities — the most likely bitstrings."""
    p = np.asarray(probs, dtype=float).ravel()
    order = np.argsort(p)[::-1][:top]
    lines = []
    for i in order:
        if p[i] < 1e-9:
            continue
        bar = "█" * max(1, int(round(p[i] * width)))
        lines.append(f"  |{format(int(i), f'0{n_qubits}b')}>  {p[i]:.4f}  {bar}")
    histogram = "\n".join(lines)
    plain = not _stream_handles_unicode() if ascii is None else ascii
    return _downgrade(histogram) if plain else histogram

qmlkit.generative

generative

Generative models — learning a distribution rather than a mapping.

Two families, and the split matters:

• Born machines (QCBM, qGAN) are implicit. Measuring the circuit samples p(x) = |<x|psi>|^2 directly, so sampling is free and scoring is not: you cannot ask such a model for p(x) of an arbitrary x without estimating it. Training therefore uses a sample-based loss — MMD, or a discriminator.
• Energy models (QBM, quantum Hopfield) are explicit. They define p(x) ∝ exp(-E(x)), so scoring is easy and sampling is hard, because the partition function sums over 2^n states.

That is the whole taxonomy, and it decides which loss you can even write down.

QCBM

QCBM(
    n_qubits: int,
    ansatz: Ansatz | None = None,
    n_layers: int = 3,
    backend: BackendLike = None,
    shots: int | None = None,
    seed: int | None = None,
)

Quantum circuit Born machine.

The circuit is the distribution: measuring it samples |<x|psi>|^2. There is no likelihood to maximise, so training minimises MMD between its samples and the data — a distance you can compute from samples alone.

Source code in src/qmlkit/generative.py

def __init__(
    self,
    n_qubits: int,
    ansatz: Ansatz | None = None,
    n_layers: int = 3,
    backend: BackendLike = None,
    shots: int | None = None,
    seed: int | None = None,
) -> None:
    self.n_qubits = n_qubits
    self.ansatz = ansatz or hardware_efficient(
        n_qubits, n_layers, entangler="cx", pattern="ring"
    )
    self.backend = backend
    self.shots = shots
    self.params_ = self.ansatz.init("uniform", seed)
    self.history_: list[float] = []

sample

sample(
    n_samples: int = 512,
    params: ArrayLike | None = None,
    seed: int | None = None,
) -> NDArray[Any]

Draw bitstrings as a (n_samples, n_qubits) array of 0/1.

Source code in src/qmlkit/generative.py

def sample(
    self, n_samples: int = 512, params: ArrayLike | None = None, seed: int | None = None
) -> npt.NDArray[Any]:
    """Draw bitstrings as a ``(n_samples, n_qubits)`` array of 0/1."""
    counts = run_counts(self.circuit(params), shots=n_samples, backend=self.backend, seed=seed)
    rows = [[int(b) for b in bits] for bits, n in counts.items() for _ in range(n)]
    return np.array(rows, dtype=int)

fit

fit(
    data: NDArray[Any],
    n_iterations: int = 100,
    gamma: float | Sequence[float] = (0.25, 1.0, 4.0),
    n_samples: int = 512,
    seed: int | None = None,
    callback: Callable[[int, NDArray[Any], float], None]
    | None = None,
) -> QCBM

Train by minimising MMD against data, using SPSA.

SPSA because the loss is a sample statistic: two circuit evaluations per step whatever the parameter count, and it tolerates the sampling noise.

Source code in src/qmlkit/generative.py

def fit(
    self,
    data: npt.NDArray[Any],
    n_iterations: int = 100,
    gamma: float | Sequence[float] = (0.25, 1.0, 4.0),
    n_samples: int = 512,
    seed: int | None = None,
    callback: Callable[[int, npt.NDArray[Any], float], None] | None = None,
) -> QCBM:
    """Train by minimising MMD against ``data``, using SPSA.

    SPSA because the loss is a sample statistic: two circuit evaluations per step
    whatever the parameter count, and it tolerates the sampling noise.
    """
    from qmlkit.gradients.spsa import minimize_spsa

    target = np.atleast_2d(np.asarray(data, dtype=float))
    rng = np.random.default_rng(seed)

    def loss(p: npt.NDArray[Any]) -> float:
        drawn = self.sample(n_samples, p, seed=int(rng.integers(1 << 31)))
        return mmd_squared(drawn, target, gamma)

    best, history = minimize_spsa(
        loss, self.params_, n_iterations=n_iterations, seed=seed, callback=callback
    )
    self.params_ = best
    self.history_ = history
    return self

score

score(
    data: NDArray[Any],
    gamma: float | Sequence[float] = (0.25, 1.0, 4.0),
    n_samples: int = 1024,
    seed: int | None = None,
) -> float

MMD² against the data — zero means the distributions match.

This is a sampled estimate, so pass seed if you need it reproducible; without one it draws on the backend's own RNG. For a deterministic comparison on a simulator use :meth:exact_distance, which needs no samples at all.

Source code in src/qmlkit/generative.py

def score(
    self,
    data: npt.NDArray[Any],
    gamma: float | Sequence[float] = (0.25, 1.0, 4.0),
    n_samples: int = 1024,
    seed: int | None = None,
) -> float:
    """MMD² against the data — zero means the distributions match.

    This is a *sampled* estimate, so pass ``seed`` if you need it reproducible;
    without one it draws on the backend's own RNG. For a deterministic comparison
    on a simulator use :meth:`exact_distance`, which needs no samples at all.
    """
    return mmd_squared(self.sample(n_samples, seed=seed), np.atleast_2d(data), gamma)

exact_distance

exact_distance(
    data: NDArray[Any],
    metric: str = "tv",
    params: Sequence[float] | None = None,
) -> float

Distance to the target distribution with no sampling.

On a simulator the model's distribution is available exactly, so progress can be measured without shot noise — which is what makes a before/after comparison trustworthy rather than a coin flip.

Source code in src/qmlkit/generative.py

def exact_distance(
    self, data: npt.NDArray[Any], metric: str = "tv", params: Sequence[float] | None = None
) -> float:
    """Distance to the target distribution with **no sampling**.

    On a simulator the model's distribution is available exactly, so progress can
    be measured without shot noise — which is what makes a before/after comparison
    trustworthy rather than a coin flip.
    """
    rows = np.atleast_2d(np.asarray(data, dtype=int))
    target = np.zeros(2**self.n_qubits)
    for row in rows:
        target[int("".join(map(str, row)), 2)] += 1
    target /= target.sum()

    model = self.probabilities(params)
    if metric == "tv":
        return total_variation(model, target)
    if metric == "kl":
        return kl_divergence(target, model)
    raise unknown("metric", metric, ("tv", "kl"))

QGAN

QGAN(
    generator: QCBM,
    discriminator: Callable[[NDArray[Any]], NDArray[Any]],
    seed: int | None = None,
)

Quantum generator, classical discriminator.

The generator is a Born machine; the discriminator is any callable scoring a batch as real. They are trained against each other, and at equilibrium the discriminator should be at chance — which is what :meth:equilibrium_gap reports.

Source code in src/qmlkit/generative.py

def __init__(
    self,
    generator: QCBM,
    discriminator: Callable[[npt.NDArray[Any]], npt.NDArray[Any]],
    seed: int | None = None,
) -> None:
    self.generator = generator
    self.discriminator = discriminator
    self.rng = np.random.default_rng(seed)
    self.history_: list[float] = []

generator_loss

generator_loss(
    params: NDArray[Any], n_samples: int = 256
) -> float

Generator wants the discriminator to call its samples real.

Source code in src/qmlkit/generative.py

def generator_loss(self, params: npt.NDArray[Any], n_samples: int = 256) -> float:
    """Generator wants the discriminator to call its samples real."""
    fake = self.generator.sample(n_samples, params, seed=int(self.rng.integers(1 << 31)))
    scores = np.asarray(self.discriminator(fake), dtype=float)
    return float(-np.log(np.clip(scores, 1e-9, 1.0)).mean())

fit_generator

fit_generator(
    n_iterations: int = 50,
    n_samples: int = 256,
    seed: int | None = None,
) -> QGAN

Train the generator against a fixed discriminator.

Source code in src/qmlkit/generative.py

def fit_generator(
    self, n_iterations: int = 50, n_samples: int = 256, seed: int | None = None
) -> QGAN:
    """Train the generator against a fixed discriminator."""
    from qmlkit.gradients.spsa import minimize_spsa

    best, history = minimize_spsa(
        lambda p: self.generator_loss(p, n_samples),
        self.generator.params_,
        n_iterations=n_iterations,
        seed=seed,
    )
    self.generator.params_ = best
    self.history_ = history
    return self

equilibrium_gap

equilibrium_gap(
    real: NDArray[Any], n_samples: int = 256
) -> float

|accuracy - 0.5| — zero when the discriminator is guessing.

Source code in src/qmlkit/generative.py

def equilibrium_gap(self, real: npt.NDArray[Any], n_samples: int = 256) -> float:
    """``|accuracy - 0.5|`` — zero when the discriminator is guessing."""
    fake = self.generator.sample(n_samples)
    real_scores = np.asarray(self.discriminator(np.atleast_2d(real)), dtype=float)
    fake_scores = np.asarray(self.discriminator(fake), dtype=float)
    accuracy = 0.5 * ((real_scores > 0.5).mean() + (fake_scores <= 0.5).mean())
    return float(abs(accuracy - 0.5))

QuantumBoltzmannMachine

QuantumBoltzmannMachine(
    n_visible: int,
    n_hidden: int = 0,
    gamma: float = 0.7,
    beta: float = 1.0,
    seed: int | None = None,
    edges: Sequence[tuple[int, int]] | None = None,
    pattern: str = "chain",
)

A transverse-field Ising model as a generative model.

The classical part is diagonal (Z fields and ZZ couplings); the transverse field Gamma * X is off-diagonal, which is what makes it quantum — and what makes the log-likelihood gradient intractable, since the model term no longer has a spins-in-number-out form. Training therefore optimises a lower bound, and grad here is the bound's clamped - model difference.

Source code in src/qmlkit/generative.py

def __init__(
    self,
    n_visible: int,
    n_hidden: int = 0,
    gamma: float = 0.7,
    beta: float = 1.0,
    seed: int | None = None,
    edges: Sequence[tuple[int, int]] | None = None,
    pattern: str = "chain",
) -> None:
    self.n_visible = n_visible
    self.n_hidden = n_hidden
    self.n_spins = n_visible + n_hidden
    self.gamma = gamma
    self.beta = beta
    rng = np.random.default_rng(seed)
    self.fields = rng.normal(0, 0.1, self.n_spins)
    # The connectivity is the model's structure, so it is an argument like any
    # other: a chain by default, "full" for an all-to-all machine, or an explicit
    # edge list for a restricted (bipartite visible/hidden) one.
    graph = (
        [(int(a), int(b)) for a, b in edges]
        if edges is not None
        else list(entangler_pairs(self.n_spins, pattern))
    )
    self.edges = graph
    self.couplings = {(a, b): float(rng.normal(0, 0.1)) for a, b in graph}

energies

energies() -> NDArray[Any]

Diagonal energy of every spin configuration.

Source code in src/qmlkit/generative.py

def energies(self) -> npt.NDArray[Any]:
    """Diagonal energy of every spin configuration."""
    configs = list(itertools.product([1, -1], repeat=self.n_spins))
    return np.array([ising_energy(c, self.fields, self.couplings) for c in configs])

probabilities

probabilities() -> NDArray[Any]

Boltzmann distribution over configurations (diagonal part only).

Source code in src/qmlkit/generative.py

def probabilities(self) -> npt.NDArray[Any]:
    """Boltzmann distribution over configurations (diagonal part only)."""
    return boltzmann(self.energies(), self.beta)[0]

visible_marginal

visible_marginal() -> NDArray[Any]

Marginalise the hidden spins away.

Source code in src/qmlkit/generative.py

def visible_marginal(self) -> npt.NDArray[Any]:
    """Marginalise the hidden spins away."""
    p = self.probabilities()
    if self.n_hidden == 0:
        return p
    return p.reshape(2**self.n_visible, 2**self.n_hidden).sum(axis=1)

grad staticmethod

grad(
    clamped: NDArray[Any], model: NDArray[Any]
) -> NDArray[Any]

data - model — the sculptor move behind every Boltzmann update.

Source code in src/qmlkit/generative.py

@staticmethod
def grad(clamped: npt.NDArray[Any], model: npt.NDArray[Any]) -> npt.NDArray[Any]:
    """``data - model`` — the sculptor move behind every Boltzmann update."""
    return np.asarray(clamped, dtype=float) - np.asarray(model, dtype=float)

free_energy

free_energy(
    entropy: float, temperature: float = 1.0
) -> float

<H> - T S — what a thermal state actually minimises.

Source code in src/qmlkit/generative.py

def free_energy(self, entropy: float, temperature: float = 1.0) -> float:
    """``<H> - T S`` — what a thermal state actually minimises."""
    return float(np.dot(self.probabilities(), self.energies()) - temperature * entropy)

QuantumHopfield

QuantumHopfield()

Associative memory: store patterns, recall the nearest by state overlap.

Recall is a fidelity comparison against each stored pattern — the same overlap a swap test estimates, which is why this belongs beside the kernel methods.

Source code in src/qmlkit/generative.py

def __init__(self) -> None:
    self.patterns_: dict[Any, npt.NDArray[Any]] = {}

store

store(
    patterns: dict[Any, Sequence[float]],
) -> QuantumHopfield

Store named patterns, normalised to unit vectors.

Source code in src/qmlkit/generative.py

def store(self, patterns: dict[Any, Sequence[float]]) -> QuantumHopfield:
    """Store named patterns, normalised to unit vectors."""
    for name, p in patterns.items():
        vec = np.asarray(p, dtype=float).ravel()
        norm = np.linalg.norm(vec)
        if norm < 1e-12:
            raise ValueError(f"pattern {name!r} is the zero vector")
        self.patterns_[name] = vec / norm
    return self

overlaps

overlaps(cue: Sequence[float]) -> dict[Any, float]

|<pattern|cue>|^2 for every stored pattern.

Source code in src/qmlkit/generative.py

def overlaps(self, cue: Sequence[float]) -> dict[Any, float]:
    """``|<pattern|cue>|^2`` for every stored pattern."""
    if not self.patterns_:
        raise ValueError("no patterns stored")
    c = np.asarray(cue, dtype=float).ravel()
    norm = np.linalg.norm(c)
    if norm < 1e-12:
        raise ValueError("cue is the zero vector")
    c = c / norm
    return {k: float(abs(np.dot(v, c)) ** 2) for k, v in self.patterns_.items()}

recall

recall(cue: Sequence[float]) -> Any

The stored pattern the cue overlaps most.

Source code in src/qmlkit/generative.py

def recall(self, cue: Sequence[float]) -> Any:
    """The stored pattern the cue overlaps most."""
    scores = self.overlaps(cue)
    return max(scores, key=lambda k: scores[k])

swap_probability staticmethod

swap_probability(overlap: float) -> float

P(anc=0) = (1 + overlap) / 2 — the swap-test readout.

Source code in src/qmlkit/generative.py

@staticmethod
def swap_probability(overlap: float) -> float:
    """``P(anc=0) = (1 + overlap) / 2`` — the swap-test readout."""
    return 0.5 + 0.5 * float(overlap)

overlap_from_probability staticmethod

overlap_from_probability(p0: float) -> float

Invert the swap-test readout.

Source code in src/qmlkit/generative.py

@staticmethod
def overlap_from_probability(p0: float) -> float:
    """Invert the swap-test readout."""
    return 2.0 * float(p0) - 1.0

gaussian_kernel

gaussian_kernel(
    a: NDArray[Any], b: NDArray[Any], gamma: float = 1.0
) -> NDArray[Any]

exp(-gamma |a - b|^2) over all pairs.

Source code in src/qmlkit/generative.py

def gaussian_kernel(
    a: npt.NDArray[Any], b: npt.NDArray[Any], gamma: float = 1.0
) -> npt.NDArray[Any]:
    """``exp(-gamma |a - b|^2)`` over all pairs."""
    x = np.atleast_2d(np.asarray(a, dtype=float))
    y = np.atleast_2d(np.asarray(b, dtype=float))
    if x.shape[1] != y.shape[1]:
        raise ValueError(f"sample widths differ: {x.shape[1]} vs {y.shape[1]}")
    sq = ((x[:, None, :] - y[None, :, :]) ** 2).sum(-1)
    return np.exp(-gamma * sq)

mmd_squared

mmd_squared(
    x: NDArray[Any],
    y: NDArray[Any],
    gamma: float | Sequence[float] = 1.0,
) -> float

Maximum mean discrepancy between two sample sets.

Zero exactly when the distributions match. A sample statistic — it never needs p(x), which is why an implicit model can be trained on it at all. Passing several gamma values averages kernels of different widths, which stops the loss going blind at one scale.

Source code in src/qmlkit/generative.py

def mmd_squared(
    x: npt.NDArray[Any], y: npt.NDArray[Any], gamma: float | Sequence[float] = 1.0
) -> float:
    """Maximum mean discrepancy between two sample sets.

    Zero exactly when the distributions match. A *sample* statistic — it never needs
    ``p(x)``, which is why an implicit model can be trained on it at all. Passing
    several ``gamma`` values averages kernels of different widths, which stops the
    loss going blind at one scale.
    """
    gammas = [gamma] if isinstance(gamma, (int, float)) else list(gamma)
    total = 0.0
    for g in gammas:
        total += float(
            gaussian_kernel(x, x, g).mean()
            + gaussian_kernel(y, y, g).mean()
            - 2 * gaussian_kernel(x, y, g).mean()
        )
    return total / len(gammas)

kl_divergence

kl_divergence(
    p: NDArray[Any], q: NDArray[Any], eps: float = 1e-12
) -> float

KL(p || q) over two discrete distributions.

Source code in src/qmlkit/generative.py

def kl_divergence(p: npt.NDArray[Any], q: npt.NDArray[Any], eps: float = 1e-12) -> float:
    """``KL(p || q)`` over two discrete distributions."""
    a = np.clip(np.asarray(p, dtype=float), eps, None)
    b = np.clip(np.asarray(q, dtype=float), eps, None)
    a, b = a / a.sum(), b / b.sum()
    return float(np.sum(a * np.log(a / b)))

total_variation

total_variation(p: NDArray[Any], q: NDArray[Any]) -> float

0.5 * sum |p - q| — bounded in [0, 1], unlike KL.

Source code in src/qmlkit/generative.py

def total_variation(p: npt.NDArray[Any], q: npt.NDArray[Any]) -> float:
    """``0.5 * sum |p - q|`` — bounded in ``[0, 1]``, unlike KL."""
    a = np.asarray(p, dtype=float)
    b = np.asarray(q, dtype=float)
    return float(0.5 * np.abs(a / a.sum() - b / b.sum()).sum())

boltzmann

boltzmann(
    energies: NDArray[Any], beta: float = 1.0
) -> tuple[NDArray[Any], float]

(p, Z) for p(x) = exp(-beta E(x)) / Z.

Source code in src/qmlkit/generative.py

def boltzmann(energies: npt.NDArray[Any], beta: float = 1.0) -> tuple[npt.NDArray[Any], float]:
    """``(p, Z)`` for ``p(x) = exp(-beta E(x)) / Z``."""
    e = np.asarray(energies, dtype=float)
    weights = np.exp(-beta * (e - e.min()))  # shift for numerical stability
    z = float(weights.sum())
    return weights / z, z

partition_function

partition_function(
    energies: NDArray[Any], beta: float = 1.0
) -> float

Z = sum exp(-beta E) — the sum over 2^n states that makes sampling hard.

Source code in src/qmlkit/generative.py

def partition_function(energies: npt.NDArray[Any], beta: float = 1.0) -> float:
    """``Z = sum exp(-beta E)`` — the sum over ``2^n`` states that makes sampling hard."""
    return float(np.exp(-beta * np.asarray(energies, dtype=float)).sum())

ising_energy

ising_energy(
    spins: Sequence[int],
    fields: NDArray[Any],
    couplings: dict[tuple[int, int], float],
) -> float

-sum b_i s_i - sum w_ij s_i s_j for spins in {+1, -1}.

Source code in src/qmlkit/generative.py

def ising_energy(
    spins: Sequence[int],
    fields: npt.NDArray[Any],
    couplings: dict[tuple[int, int], float],
) -> float:
    """``-sum b_i s_i - sum w_ij s_i s_j`` for spins in ``{+1, -1}``."""
    s = np.asarray(spins, dtype=float)
    energy = -float(np.dot(np.asarray(fields, dtype=float), s))
    for (i, j), w in couplings.items():
        energy -= float(w) * float(s[i]) * float(s[j])
    return energy

qmlkit.shadows

Classical shadows: many observables from few measurements.

shadows

Classical shadows — estimate many observables from few measurements.

Huang, Kueng & Preskill (2020). Measure in a randomly chosen basis each shot, invert the resulting depolarising channel, and the collection of snapshots predicts :math:M observables to fixed accuracy from :math:O(\log M) measurements — rather than measuring each one separately.

On an exact simulator this buys nothing: shots=None returns every observable exactly. It earns its place because measurement cost is what binds on hardware, and because it makes that cost visible: :func:shadow_shot_cost against :func:qmlkit.kernels.kernel_shot_cost is the comparison worth looking at before committing to a device run.

shadow = ClassicalShadow(spec, n_snapshots=2000, seed=0)
shadow.expectation(qk.Z(0) + 0.5 * qk.ZZ(0, 2))

ClassicalShadow

ClassicalShadow(
    spec: CircuitSpec,
    n_snapshots: int = 1000,
    seed: int | None = None,
    backend: BackendLike = None,
)

A set of randomised single-qubit measurements, and what they predict.

Each snapshot picks a random Pauli basis per qubit, measures once, and stores (basis, outcome). An observable's estimate then averages over only the snapshots whose bases happen to match its support — which is why the cost grows with the observable's locality, not with how many observables you ask for.

Source code in src/qmlkit/shadows.py

def __init__(
    self,
    spec: CircuitSpec,
    n_snapshots: int = 1000,
    seed: int | None = None,
    backend: BackendLike = None,
) -> None:
    if not spec.is_bound:
        raise ValueError("bind the circuit before taking shadows")
    self.spec = spec
    self.n_qubits = spec.n_qubits
    self.n_snapshots = n_snapshots
    self.backend = backend
    rng = np.random.default_rng(seed)
    #: (n_snapshots, n_qubits) of basis indices, and of +-1 outcomes
    self.bases = rng.integers(0, 3, size=(n_snapshots, self.n_qubits))
    self.outcomes = np.empty((n_snapshots, self.n_qubits), dtype=np.int8)
    self._collect(rng)

expectation

expectation(obs: Observable) -> float

Estimate <O> from the stored snapshots.

Source code in src/qmlkit/shadows.py

def expectation(self, obs: Observable) -> float:
    """Estimate ``<O>`` from the stored snapshots."""
    total = 0.0
    for term in as_sum(obs).terms:
        total += float(term.coeff.real) * self._term_estimate(term)
    return total

shadow_shot_cost

shadow_shot_cost(
    locality: int, n_observables: int, epsilon: float = 0.1
) -> int

Snapshots for :math:M observables of given locality to accuracy epsilon.

The headline scaling: :math:O(3^k \log M / \epsilon^2) — logarithmic in how many observables you want, exponential only in their locality. Measuring each one separately is instead linear in M, which is the trade this whole method makes.

Source code in src/qmlkit/shadows.py

def shadow_shot_cost(locality: int, n_observables: int, epsilon: float = 0.1) -> int:
    r"""Snapshots for :math:`M` observables of given locality to accuracy ``epsilon``.

    The headline scaling: :math:`O(3^k \log M / \epsilon^2)` — **logarithmic** in how
    many observables you want, exponential only in their locality. Measuring each one
    separately is instead linear in ``M``, which is the trade this whole method makes.
    """
    if locality < 1 or n_observables < 1:
        raise ValueError("locality and n_observables must be at least 1")
    return int(np.ceil(3**locality * np.log(2 * n_observables) / epsilon**2))

qmlkit.utils.shots

shots

Shot-budget arithmetic.

Simulator-only means shots are opt-in, not mandatory. When they are on, every sampled number should be reportable with its uncertainty — that is what makes "would this survive on a real device?" an answerable question rather than a guess.

variance

variance(z: float) -> float

Single-shot variance of a +-1 observable with mean z: 1 - z**2.

Source code in src/qmlkit/utils/shots.py

def variance(z: float) -> float:
    """Single-shot variance of a +-1 observable with mean ``z``: ``1 - z**2``."""
    return float(1.0 - np.clip(z, -1.0, 1.0) ** 2)

standard_error

standard_error(
    z: float, shots: int, scale: float = 1.0
) -> float

Standard error of a single Pauli term's expectation over shots samples.

scale is the term's coefficient: (cP)^2 = c^2 I, so the variance is c^2 - z^2 and the error scales with |c|.

This formula is only correct for one term. A sum needs <O^2>, which is a different measurement and not recoverable from <O> — feeding a sum in here gives an error bar that is too tight, too loose, or exactly zero once |<O>| reaches |c|. :func:~qmlkit.core.execute.expectation checks the observable before calling this.

Source code in src/qmlkit/utils/shots.py

def standard_error(z: float, shots: int, scale: float = 1.0) -> float:
    """Standard error of a **single** Pauli term's expectation over ``shots`` samples.

    ``scale`` is the term's coefficient: ``(cP)^2 = c^2 I``, so the variance is
    ``c^2 - z^2`` and the error scales with ``|c|``.

    This formula is only correct for one term. A sum needs ``<O^2>``, which is a
    different measurement and not recoverable from ``<O>`` — feeding a sum in here
    gives an error bar that is too tight, too loose, or exactly zero once ``|<O>|``
    reaches ``|c|``. :func:`~qmlkit.core.execute.expectation` checks the observable
    before calling this.
    """
    if shots <= 0:
        raise ValueError("shots must be positive")
    magnitude = abs(float(scale))
    spread = max(magnitude**2 - float(z) ** 2, 0.0)
    return float(np.sqrt(spread / shots))

shots_for_precision

shots_for_precision(eps: float, z: float = 0.0) -> int

Shots needed to reach standard error eps — the 1/eps**2 price.

Source code in src/qmlkit/utils/shots.py

def shots_for_precision(eps: float, z: float = 0.0) -> int:
    """Shots needed to reach standard error ``eps`` — the ``1/eps**2`` price."""
    if eps <= 0:
        raise ValueError("eps must be positive")
    return int(np.ceil(variance(z) / eps**2))

p0_from_z

p0_from_z(z: float) -> float

P(0) from .

Source code in src/qmlkit/utils/shots.py

def p0_from_z(z: float) -> float:
    """P(0) from <Z>."""
    return float((1.0 + z) / 2.0)

z_from_p0

z_from_p0(p0: float) -> float

from P(0).

Source code in src/qmlkit/utils/shots.py

def z_from_p0(p0: float) -> float:
    """<Z> from P(0)."""
    return float(2.0 * p0 - 1.0)

runtime_estimate

runtime_estimate(shots: int, rate_hz: float) -> float

Wall-clock seconds for a shot budget at a given sampling rate.

Source code in src/qmlkit/utils/shots.py

def runtime_estimate(shots: int, rate_hz: float) -> float:
    """Wall-clock seconds for a shot budget at a given sampling rate."""
    if rate_hz <= 0:
        raise ValueError("rate_hz must be positive")
    return float(shots / rate_hz)