Skip to content

Kernels

Three overlap estimators, Gram matrices that stay positive semi-definite, and the models on top.

qmlkit.kernels.estimators

estimators

The three ways to read a kernel off a circuit.

A quantum kernel is an overlap: :math:k(x, x') = |\langle \phi(x')|\phi(x)\rangle|^2. Three estimators get at it, and they are not interchangeable:

============ ========== ==================== ========================================= Estimator Qubits Depth Gives you ============ ========== ==================== ========================================= inversion n 2 x feature map |<.>|^2 — the magnitude swap test 2n + 1 1 map + n CSWAPs |<.>|^2, from an ancilla Hadamard n + 1 controlled map Re<.> — signed ============ ========== ==================== =========================================

The inversion (compute-uncompute) test is the default: fewest qubits, no ancilla, no controlled gates. The swap test earns its extra register when you already hold two states and cannot rebuild one. The Hadamard test is the only one that keeps the sign of the inner product — magnitude estimators map +1/2 and -1/2 to the same number.

inversion_circuit 

inversion_circuit(
    fmap: FeatureMap,
    x: Sequence[float],
    xp: Sequence[float],
) -> CircuitSpec

U(x) followed by U(x')†. P(all zeros) is the kernel.

Source code in src/qmlkit/kernels/estimators.py 
43
44
45
def inversion_circuit(fmap: FeatureMap, x: Sequence[float], xp: Sequence[float]) -> CircuitSpec:
    """``U(x)`` followed by ``U(x')†``. P(all zeros) *is* the kernel."""
    return fmap.build(x).compose(fmap.adjoint(xp), param_offset=0)

fidelity_kernel 

fidelity_kernel(
    fmap: FeatureMap,
    x: Sequence[float],
    xp: Sequence[float],
    shots: int | None = None,
    backend: BackendLike = None,
    seed: int | None = None,
) -> float

k(x, x') by the compute-uncompute test — the default estimator.

With shots=None this reads the exact all-zeros probability. With a shot budget it counts the all-zeros outcomes, which is what a device would do.

Source code in src/qmlkit/kernels/estimators.py 
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
def fidelity_kernel(
    fmap: FeatureMap,
    x: Sequence[float],
    xp: Sequence[float],
    shots: int | None = None,
    backend: BackendLike = None,
    seed: int | None = None,
) -> float:
    """``k(x, x')`` by the compute-uncompute test — the default estimator.

    With ``shots=None`` this reads the exact all-zeros probability. With a shot
    budget it counts the all-zeros outcomes, which is what a device would do.
    """
    spec = inversion_circuit(fmap, x, xp)
    if shots is None:
        return float(probabilities(spec, backend=backend)[0])
    counts = run_counts(spec, shots=shots, backend=backend, seed=seed)
    zeros = "0" * spec.n_qubits
    return counts.get(zeros, 0) / shots

swap_readout 

swap_readout(p_ancilla_zero: float) -> float

Invert the swap-test readout: k = 2 P(anc=0) - 1.

Source code in src/qmlkit/kernels/estimators.py 
69
70
71
def swap_readout(p_ancilla_zero: float) -> float:
    """Invert the swap-test readout: ``k = 2 P(anc=0) - 1``."""
    return 2.0 * float(p_ancilla_zero) - 1.0

swap_probability 

swap_probability(k: float) -> float

Forward direction: P(anc=0) = (1 + k) / 2. Orthogonal states give a fair coin.

Source code in src/qmlkit/kernels/estimators.py 
74
75
76
def swap_probability(k: float) -> float:
    """Forward direction: ``P(anc=0) = (1 + k) / 2``. Orthogonal states give a fair coin."""
    return (1.0 + float(k)) / 2.0

swap_test_kernel 

swap_test_kernel(
    fmap: FeatureMap,
    x: Sequence[float],
    xp: Sequence[float],
    shots: int | None = None,
    backend: BackendLike = None,
    seed: int | None = None,
) -> float

k(x, x') by the swap test — two registers plus one ancilla.

Costs 2n + 1 qubits against the inversion test's n, and needs a CSWAP per qubit pair. Worth it only when you genuinely hold two states already.

Source code in src/qmlkit/kernels/estimators.py 
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
def swap_test_kernel(
    fmap: FeatureMap,
    x: Sequence[float],
    xp: Sequence[float],
    shots: int | None = None,
    backend: BackendLike = None,
    seed: int | None = None,
) -> float:
    """``k(x, x')`` by the swap test — two registers plus one ancilla.

    Costs ``2n + 1`` qubits against the inversion test's ``n``, and needs a CSWAP
    per qubit pair. Worth it only when you genuinely hold two states already.
    """
    a, b = fmap.build(x), fmap.build(xp)
    n = fmap.n_qubits
    ancilla = 2 * n

    ops: list[Op] = []
    ops += [Op(op.gate, tuple(q for q in op.qubits), op.params) for op in a.ops]
    ops += [Op(op.gate, tuple(q + n for q in op.qubits), op.params) for op in b.ops]
    ops.append(Op("h", (ancilla,)))
    for i in range(n):
        ops += _controlled_swap(ancilla, i, i + n)
    ops.append(Op("h", (ancilla,)))

    spec = CircuitSpec(2 * n + 1, tuple(ops), 0)
    if shots is None:
        probs = probabilities(spec, backend=backend)
        # ancilla is the least significant bit (qubit 0 is most significant)
        p0 = float(probs[::2].sum())
    else:
        counts = run_counts(spec, shots=shots, backend=backend, seed=seed)
        p0 = sum(v for k, v in counts.items() if k[ancilla] == "0") / shots
    return swap_readout(p0)

hadamard_test 

hadamard_test(
    fmap: FeatureMap,
    x: Sequence[float],
    xp: Sequence[float],
    part: str = "real",
    shots: int | None = None,
    backend: BackendLike = None,
    seed: int | None = None,
) -> float

Re<phi(x')|phi(x)> (or the imaginary part) — the signed inner product.

The only estimator that distinguishes +1/2 from -1/2; the magnitude ones map both to 1/4. The price is an ancilla and a controlled feature map, which roughly doubles depth and wants all-to-all connectivity to the ancilla — which is why it is rarely the right choice on hardware.

Source code in src/qmlkit/kernels/estimators.py 
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
def hadamard_test(
    fmap: FeatureMap,
    x: Sequence[float],
    xp: Sequence[float],
    part: str = "real",
    shots: int | None = None,
    backend: BackendLike = None,
    seed: int | None = None,
) -> float:
    """``Re<phi(x')|phi(x)>`` (or the imaginary part) — the **signed** inner product.

    The only estimator that distinguishes ``+1/2`` from ``-1/2``; the magnitude ones
    map both to ``1/4``. The price is an ancilla and a *controlled* feature map,
    which roughly doubles depth and wants all-to-all connectivity to the ancilla —
    which is why it is rarely the right choice on hardware.
    """
    if part not in ("real", "imag"):
        raise unknown("part", part, ("real", "imag"))
    n = fmap.n_qubits
    ancilla = n

    qc = QCircuit(n + 1)
    qc.h(ancilla)
    if part == "imag":
        qc.sdg(ancilla)
    spec = qc.to_spec()

    body = fmap.build(x).compose(fmap.adjoint(xp), param_offset=0)
    ops = list(spec.ops)
    for op in body.ops:
        ops.extend(_controlled(op, ancilla))
    ops.append(Op("h", (ancilla,)))

    full = CircuitSpec(n + 1, tuple(ops), 0)
    if shots is None:
        probs = probabilities(full, backend=backend)
        p0 = float(probs[::2].sum())
    else:
        counts = run_counts(full, shots=shots, backend=backend, seed=seed)
        p0 = sum(v for k, v in counts.items() if k[ancilla] == "0") / shots
    return 2.0 * p0 - 1.0

qmlkit.kernels.matrix

matrix

Gram matrices, PSD repair, and the diagnostics that say whether any of it will work.

Filling a Gram matrix is the expensive half of a quantum kernel method: m(m-1)/2 circuit evaluations for a training set of size m, since the diagonal is exactly 1 and the matrix is symmetric. :func:kernel_matrix exploits both.

Shot noise breaks positive semi-definiteness. Every entry is an estimate, so the estimated Gram matrix can have small negative eigenvalues even though the true one cannot — and an SVM solver will either refuse it or return nonsense. The repair functions here project back onto the PSD cone.

Exponential concentration is the real limit. As the feature map widens, distinct inputs produce states whose overlaps all collapse toward the same value, at a rate around 2^-n. Resolving that against shot noise costs about 4^n shots. The diagnostics report both, because a kernel method that has concentrated looks like a model that simply does not learn.

QuantumKernel 

QuantumKernel(
    feature_map: FeatureMap,
    estimator: str = "inversion",
    shots: int | None = None,
    backend: BackendLike = None,
    bandwidth: float = 1.0,
    seed: int | None = None,
    cache: bool = True,
)

A feature map, as a kernel you can hand to any kernel method.

kernel = QuantumKernel(qk.ZZFeatureMap(2)) K = kernel(X) # training Gram matrix K_test = kernel(X_test, X) # rectangular, test against train

Arguments

estimator How the overlap is measured. "inversion" (the default) runs the compute-uncompute circuit and reads the all-zeros probability; "swap" uses a swap test; "hadamard" runs two Hadamard tests and adds the squares of the real and imaginary parts, because one Hadamard test measures one component of a complex overlap and the kernel is its modulus. All three agree on a simulator — :attr:n_evaluations is what differs, and on a device so do the width and the connectivity each one needs. shots None reads the exact probability. A budget samples it, which is what a device does — and a sampled kernel is not positive semi-definite by construction, so pair it with :func:threshold_matrix. bandwidth The first thing to try when a kernel has concentrated. Every feature vector is scaled by this before encoding, so it sets how far apart two points are in the feature map rather than in the data. At the default 1.0 a fidelity kernel over a wide register drives every off-diagonal entry toward the same small number — every pair of points looks equally dissimilar, the Gram matrix approaches the identity, and no amount of training recovers what the encoding threw away. Shrinking the bandwidth (0.1-0.5 is the usual range) compresses the data into a smaller region of state space and pulls the off-diagonals back apart. :func:concentration_report measures whether you have the problem, and qk.diagnose(K) names it as KERNEL_CONCENTRATED. The alternative fix is a projected kernel, which survives width by measuring local reduced states instead — see the kernels tutorial for when each applies. cache Memoises pair evaluations, which matters because a Gram matrix asks for the same circuit many times. n_evaluations counts the circuits actually run.

Source code in src/qmlkit/kernels/matrix.py 
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
def __init__(
    self,
    feature_map: FeatureMap,
    estimator: str = "inversion",
    shots: int | None = None,
    backend: BackendLike = None,
    bandwidth: float = 1.0,
    seed: int | None = None,
    cache: bool = True,
) -> None:
    self.feature_map = feature_map
    self.estimator = estimator
    self.shots = shots
    self.backend = backend
    self.bandwidth = bandwidth
    self.seed = seed
    self.cache = cache
    self._cache: dict[tuple[Any, ...], float] = {}
    self._evaluations = 0
n_evaluations property 
n_evaluations: int

Circuits actually run — cache hits do not count.

evaluate 
evaluate(x: Sequence[float], xp: Sequence[float]) -> float

One kernel entry, with the bandwidth rescaling applied.

Source code in src/qmlkit/kernels/matrix.py 
222
223
224
225
226
227
228
229
230
231
232
233
def evaluate(self, x: Sequence[float], xp: Sequence[float]) -> float:
    """One kernel entry, with the bandwidth rescaling applied."""
    a = self.bandwidth * np.asarray(x, dtype=float)
    b = self.bandwidth * np.asarray(xp, dtype=float)
    if not self.cache:
        return self._estimate(a, b)
    # the kernel is symmetric, so k(a, b) and k(b, a) share one cache entry --
    # which halves the evaluations a rectangular test matrix needs
    key = self._cache_key(a, b)
    if key not in self._cache:
        self._cache[key] = self._estimate(a, b)
    return self._cache[key]

square_kernel_matrix 

square_kernel_matrix(
    X: NDArray[Any],
    kernel: KernelFn,
    assume_unit_diagonal: bool = True,
) -> NDArray[Any]

Symmetric Gram matrix, evaluating only the upper triangle.

m(m-1)/2 evaluations instead of m^2. assume_unit_diagonal sets k(x, x) = 1 without measuring it, which is exact for a fidelity kernel and saves m more evaluations.

Source code in src/qmlkit/kernels/matrix.py 
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
def square_kernel_matrix(
    X: npt.NDArray[Any], kernel: KernelFn, assume_unit_diagonal: bool = True
) -> npt.NDArray[Any]:
    """Symmetric Gram matrix, evaluating only the upper triangle.

    ``m(m-1)/2`` evaluations instead of ``m^2``. ``assume_unit_diagonal`` sets
    ``k(x, x) = 1`` without measuring it, which is exact for a fidelity kernel and
    saves ``m`` more evaluations.
    """
    rows = np.atleast_2d(np.asarray(X, dtype=float))
    m = rows.shape[0]
    out = np.eye(m) if assume_unit_diagonal else np.zeros((m, m))
    if not assume_unit_diagonal:
        for i in range(m):
            out[i, i] = kernel(rows[i], rows[i])
    # this is the path that takes hours on a few hundred points, so it is the one
    # worth being able to watch
    with progress_task("kernel gram", m * (m - 1) // 2) as tracked:
        for i in range(m):
            for j in range(i + 1, m):
                out[i, j] = out[j, i] = kernel(rows[i], rows[j])
                tracked.advance()
    return out

kernel_matrix 

kernel_matrix(
    X: NDArray[Any],
    Y: NDArray[Any] | None = None,
    kernel: KernelFn | None = None,
) -> NDArray[Any]

Gram matrix of X against Y (or itself, exploiting symmetry).

Source code in src/qmlkit/kernels/matrix.py 
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
def kernel_matrix(
    X: npt.NDArray[Any], Y: npt.NDArray[Any] | None = None, kernel: KernelFn | None = None
) -> npt.NDArray[Any]:
    """Gram matrix of ``X`` against ``Y`` (or itself, exploiting symmetry)."""
    if kernel is None:
        raise ValueError("kernel_matrix needs a kernel callable")
    if Y is None:
        return square_kernel_matrix(X, kernel)
    a = np.atleast_2d(np.asarray(X, dtype=float))
    b = np.atleast_2d(np.asarray(Y, dtype=float))
    with progress_task("kernel gram", a.shape[0] * b.shape[0]) as tracked:
        out = np.empty((a.shape[0], b.shape[0]), dtype=float)
        for i, u in enumerate(a):
            for j, v in enumerate(b):
                out[i, j] = kernel(u, v)
                tracked.advance()
    return out

is_psd 

is_psd(K: NDArray[Any], tol: float = 1e-09) -> bool

True if every eigenvalue is non-negative to within tol.

Source code in src/qmlkit/kernels/matrix.py 
347
348
349
def is_psd(K: npt.NDArray[Any], tol: float = 1e-9) -> bool:
    """True if every eigenvalue is non-negative to within ``tol``."""
    return min_eigenvalue(K) >= -abs(tol)

threshold_matrix 

threshold_matrix(K: NDArray[Any]) -> NDArray[Any]

Clip negative eigenvalues to zero — the standard projection onto the cone.

Source code in src/qmlkit/kernels/matrix.py 
352
353
354
355
def threshold_matrix(K: npt.NDArray[Any]) -> npt.NDArray[Any]:
    """Clip negative eigenvalues to zero — the standard projection onto the cone."""
    vals, vecs = np.linalg.eigh(np.asarray(K, dtype=float))
    return (vecs * np.clip(vals, 0.0, None)) @ vecs.T

displace_matrix 

displace_matrix(K: NDArray[Any]) -> NDArray[Any]

Shift the whole spectrum up until it is non-negative.

Keeps every eigenvector's relative weight, unlike thresholding, at the cost of inflating the diagonal.

Source code in src/qmlkit/kernels/matrix.py 
358
359
360
361
362
363
364
365
366
def displace_matrix(K: npt.NDArray[Any]) -> npt.NDArray[Any]:
    """Shift the whole spectrum up until it is non-negative.

    Keeps every eigenvector's relative weight, unlike thresholding, at the cost of
    inflating the diagonal.
    """
    arr = np.asarray(K, dtype=float)
    low = min_eigenvalue(arr)
    return arr if low >= 0 else arr + abs(low) * np.eye(arr.shape[0])

flip_matrix 

flip_matrix(K: NDArray[Any]) -> NDArray[Any]

Take the absolute value of each eigenvalue.

Source code in src/qmlkit/kernels/matrix.py 
369
370
371
372
def flip_matrix(K: npt.NDArray[Any]) -> npt.NDArray[Any]:
    """Take the absolute value of each eigenvalue."""
    vals, vecs = np.linalg.eigh(np.asarray(K, dtype=float))
    return (vecs * np.abs(vals)) @ vecs.T

closest_psd_matrix 

closest_psd_matrix(
    K: NDArray[Any], method: str = "threshold"
) -> NDArray[Any]

Nearest PSD matrix by the named method.

Source code in src/qmlkit/kernels/matrix.py 
375
376
377
378
379
380
381
def closest_psd_matrix(K: npt.NDArray[Any], method: str = "threshold") -> npt.NDArray[Any]:
    """Nearest PSD matrix by the named method."""
    fns = {"threshold": threshold_matrix, "displace": displace_matrix, "flip": flip_matrix}
    try:
        return fns[method](K)
    except KeyError:
        raise unknown("method", method, ("threshold", "displace", "flip")) from None

center_kernel 

center_kernel(K: NDArray[Any]) -> NDArray[Any]

Centre the induced feature space at the origin.

Source code in src/qmlkit/kernels/matrix.py 
384
385
386
387
388
389
def center_kernel(K: npt.NDArray[Any]) -> npt.NDArray[Any]:
    """Centre the induced feature space at the origin."""
    arr = np.asarray(K, dtype=float)
    m = arr.shape[0]
    ones = np.ones((m, m)) / m
    return arr - ones @ arr - arr @ ones + ones @ arr @ ones

normalize_kernel 

normalize_kernel(K: NDArray[Any]) -> NDArray[Any]

Rescale to a unit diagonal — the cosine of the feature-space angle.

Source code in src/qmlkit/kernels/matrix.py 
392
393
394
395
396
def normalize_kernel(K: npt.NDArray[Any]) -> npt.NDArray[Any]:
    """Rescale to a unit diagonal — the cosine of the feature-space angle."""
    arr = np.asarray(K, dtype=float)
    d = np.sqrt(np.clip(np.diag(arr), 1e-15, None))
    return arr / np.outer(d, d)

target_alignment 

target_alignment(
    K: NDArray[Any], y: NDArray[Any], rescale: bool = True
) -> float

Kernel-target alignment: how much the Gram matrix looks like the labels.

<K, yy^T>_F / (||K||_F ||yy^T||_F) in [-1, 1]. This is the objective you maximise to train a feature map, and a cheap way to compare candidates without fitting an SVM to each.

Source code in src/qmlkit/kernels/matrix.py 
402
403
404
405
406
407
408
409
410
411
412
413
414
415
def target_alignment(K: npt.NDArray[Any], y: npt.NDArray[Any], rescale: bool = True) -> float:
    """Kernel-target alignment: how much the Gram matrix looks like the labels.

    ``<K, yy^T>_F / (||K||_F ||yy^T||_F)`` in ``[-1, 1]``. This is the objective you
    maximise to *train* a feature map, and a cheap way to compare candidates without
    fitting an SVM to each.
    """
    arr = np.asarray(K, dtype=float)
    labels = np.asarray(y, dtype=float).ravel()
    if rescale and set(np.unique(labels)) <= {0.0, 1.0}:
        labels = 2 * labels - 1  # {0,1} -> {-1,+1}
    target = np.outer(labels, labels)
    denom = np.linalg.norm(arr) * np.linalg.norm(target)
    return float(np.sum(arr * target) / denom) if denom > 0 else 0.0

kernel_shot_cost 

kernel_shot_cost(
    m: int, shots: int, include_diagonal: bool = False
) -> int

Total shots to fill an m x m Gram matrix.

Source code in src/qmlkit/kernels/matrix.py 
418
419
420
421
def kernel_shot_cost(m: int, shots: int, include_diagonal: bool = False) -> int:
    """Total shots to fill an ``m x m`` Gram matrix."""
    entries = m * (m - 1) // 2 + (m if include_diagonal else 0)
    return entries * shots

kernel_spread 

kernel_spread(n_qubits: int) -> float

Rough off-diagonal spread of a concentrated kernel: 2^-n.

Source code in src/qmlkit/kernels/matrix.py 
424
425
426
def kernel_spread(n_qubits: int) -> float:
    """Rough off-diagonal spread of a concentrated kernel: ``2^-n``."""
    return float(2.0**-n_qubits)

shots_to_resolve 

shots_to_resolve(n_qubits: int) -> int

Shots needed to see a 2^-n signal above 1/sqrt(N) noise: about 4^n.

Source code in src/qmlkit/kernels/matrix.py 
429
430
431
def shots_to_resolve(n_qubits: int) -> int:
    """Shots needed to see a ``2^-n`` signal above ``1/sqrt(N)`` noise: about ``4^n``."""
    return int(4**n_qubits)

concentration_report 

concentration_report(
    K: NDArray[Any], n_qubits: int, shots: int | None = None
) -> dict[str, Any]

Is this Gram matrix telling you anything, or has it concentrated?

A concentrated kernel has near-identical off-diagonal entries: every pair of inputs looks equally similar, so no model built on it can separate them.

Source code in src/qmlkit/kernels/matrix.py 
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
def concentration_report(
    K: npt.NDArray[Any], n_qubits: int, shots: int | None = None
) -> dict[str, Any]:
    """Is this Gram matrix telling you anything, or has it concentrated?

    A concentrated kernel has near-identical off-diagonal entries: every pair of
    inputs looks equally similar, so no model built on it can separate them.
    """
    arr = np.asarray(K, dtype=float)
    off = arr[~np.eye(arr.shape[0], dtype=bool)]
    spread = float(off.std())
    noise = float(np.sqrt(0.25 / shots)) if shots else 0.0
    return {
        "off_diagonal_mean": float(off.mean()),
        "off_diagonal_std": spread,
        "predicted_spread": kernel_spread(n_qubits),
        "shot_noise": noise,
        "resolvable": bool(spread > noise) if shots else True,
        "shots_to_resolve": shots_to_resolve(n_qubits),
        "min_eigenvalue": min_eigenvalue(arr),
        "is_psd": is_psd(arr),
    }

geometric_difference 

geometric_difference(
    k_quantum: NDArray[Any], k_classical: NDArray[Any]
) -> float

g(K_C || K_Q) — the statistic that says whether quantum could help.

.. math::

g = \sqrt{\lVert \sqrt{K_Q}\, K_C^{-1} \sqrt{K_Q} \rVert_\infty}

The number to compare it against is sqrt(N), for N samples — that is the threshold in Huang et al. (2021), and it is the caller's to apply. A g well below sqrt(N) says the classical kernel already sees everything the quantum one does, so no separation is available whatever a later accuracy table claims. A g at or above it says a separation is possible, not that one exists.

Identical kernels give exactly 1. The statistic is scale-sensitive, so it carries the paper's meaning only when the two kernels are normalised alike — which any two with a unit diagonal are, a fidelity kernel and an RBF kernel included. center_kernel or normalize_kernel will put an odd one right.

Source code in src/qmlkit/kernels/matrix.py 
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
def geometric_difference(k_quantum: npt.NDArray[Any], k_classical: npt.NDArray[Any]) -> float:
    r"""``g(K_C || K_Q)`` — the statistic that says whether quantum *could* help.

    .. math::

        g = \sqrt{\lVert \sqrt{K_Q}\, K_C^{-1} \sqrt{K_Q} \rVert_\infty}

    **The number to compare it against is** ``sqrt(N)``, for ``N`` samples — that is
    the threshold in Huang et al. (2021), and it is the caller's to apply. A ``g``
    well below ``sqrt(N)`` says the classical kernel already sees everything the
    quantum one does, so no separation is available whatever a later accuracy table
    claims. A ``g`` at or above it says a separation is *possible*, not that one
    exists.

    Identical kernels give exactly ``1``. The statistic is scale-sensitive, so it
    carries the paper's meaning only when the two kernels are normalised alike —
    which any two with a unit diagonal are, a fidelity kernel and an RBF kernel
    included. ``center_kernel`` or ``normalize_kernel`` will put an odd one right.
    """
    kq = np.asarray(k_quantum, dtype=float)
    kc = np.asarray(k_classical, dtype=float)
    if kq.shape != kc.shape:
        raise ValueError(f"kernels have different shapes: {kq.shape} vs {kc.shape}")
    n = kq.shape[0]
    sqrt_kq = _sqrtm_psd(kq)
    kc_inv = np.linalg.pinv(kc + 1e-12 * np.eye(n))
    m = sqrt_kq @ kc_inv @ sqrt_kq
    return float(np.sqrt(np.linalg.norm(m, ord=2)))

qmlkit.kernels.models

models

Kernel models: sklearn estimators, a trainable kernel, and projected kernels.

The division of labour a quantum kernel method rests on: the quantum part fills the Gram matrix, and the classical part solves a convex problem on it. That means QSVC is a real SVM — same convergence guarantees, same solver — with one matrix supplied from a circuit.

QSVC 

QSVC(
    feature_map: FeatureMap, C: float = 1.0, **kwargs: Any
)

Bases: _KernelEstimator

Source code in src/qmlkit/kernels/models.py 
122
123
124
125
126
def __init__(self, feature_map: FeatureMap, C: float = 1.0, **kwargs: Any) -> None:
    super().__init__(feature_map, **kwargs)
    self.C = C  # stored under its own name, so get_params/clone can see it
    svm = _require_sklearn("QSVC")
    self._svm = svm.SVC(kernel="precomputed", C=C, **self.solver_kwargs)

QSVR 

QSVR(
    feature_map: FeatureMap,
    C: float = 1.0,
    epsilon: float = 0.1,
    **kwargs: Any,
)

Bases: _KernelEstimator

Quantum-kernel support vector regressor.

Source code in src/qmlkit/kernels/models.py 
137
138
139
140
141
142
143
144
def __init__(
    self, feature_map: FeatureMap, C: float = 1.0, epsilon: float = 0.1, **kwargs: Any
) -> None:
    super().__init__(feature_map, **kwargs)
    self.C = C
    self.epsilon = epsilon
    svm = _require_sklearn("QSVR")
    self._svm = svm.SVR(kernel="precomputed", C=C, epsilon=epsilon, **self.solver_kwargs)
score 
score(X: NDArray[Any], y: NDArray[Any]) -> float

R^2.

Source code in src/qmlkit/kernels/models.py 
146
147
148
149
150
151
152
def score(self, X: npt.NDArray[Any], y: npt.NDArray[Any]) -> float:
    """R^2."""
    pred = np.asarray(self.predict(X)).ravel()
    truth = np.asarray(y, dtype=float).ravel()
    ss_res = float(((truth - pred) ** 2).sum())
    ss_tot = float(((truth - truth.mean()) ** 2).sum())
    return 1.0 - ss_res / ss_tot if ss_tot > 0 else 0.0

NearestFidelityClassifier 

NearestFidelityClassifier(
    feature_map: FeatureMap,
    shots: int | None = None,
    backend: BackendLike = None,
)

Classify by fidelity to each class centroid — no solver, no sklearn.

The simplest quantum classifier there is: encode every training point, average within each class, and predict whichever class anchor a new point overlaps most.

Source code in src/qmlkit/kernels/models.py 
162
163
164
165
166
167
168
169
170
171
172
def __init__(
    self,
    feature_map: FeatureMap,
    shots: int | None = None,
    backend: BackendLike = None,
) -> None:
    self.feature_map = feature_map
    self.shots = shots
    self.backend = backend
    self.classes_: npt.NDArray[Any] | None = None
    self.anchors_: dict[Any, npt.NDArray[Any]] = {}

TrainableKernel 

TrainableKernel(
    feature_map_factory: Any,
    n_params: int,
    shots: int | None = None,
    backend: BackendLike = None,
)

Train the feature map itself by maximising kernel-target alignment.

A fixed feature map is a guess. Alignment gives a differentiable score for how well a kernel matches the labels, so the embedding's own parameters can be optimised before any classifier is fitted — usually a bigger win than tuning the classifier afterwards.

Source code in src/qmlkit/kernels/models.py 
216
217
218
219
220
221
222
223
224
225
226
227
228
def __init__(
    self,
    feature_map_factory: Any,
    n_params: int,
    shots: int | None = None,
    backend: BackendLike = None,
) -> None:
    self.factory = feature_map_factory
    self.n_params = n_params
    self.shots = shots
    self.backend = backend
    self.params_: npt.NDArray[Any] | None = None
    self.history_: list[float] = []
fit 
fit(
    X: NDArray[Any],
    y: NDArray[Any],
    n_iterations: int = 40,
    theta0: Sequence[float] | None = None,
    seed: int | None = None,
) -> TrainableKernel

Maximise alignment with SPSA — two evaluations per step, any parameter count.

Source code in src/qmlkit/kernels/models.py 
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
def fit(
    self,
    X: npt.NDArray[Any],
    y: npt.NDArray[Any],
    n_iterations: int = 40,
    theta0: Sequence[float] | None = None,
    seed: int | None = None,
) -> TrainableKernel:
    """Maximise alignment with SPSA — two evaluations per step, any parameter count."""
    from qmlkit.gradients.spsa import minimize_spsa

    rng = np.random.default_rng(seed)
    start = (
        np.asarray(theta0, dtype=float)
        if theta0 is not None
        else rng.uniform(0.5, 1.5, self.n_params)
    )

    def loss(p: npt.NDArray[Any]) -> float:
        return -self.alignment(p, X, y)  # minimise the negative

    # minimize_spsa returns the *final* iterate, not the best one seen, so a short
    # run can legitimately end below where it started. history_ is the alignment
    # trajectory and its last entry always corresponds to params_.
    final, history = minimize_spsa(loss, start, n_iterations=n_iterations, seed=seed)
    self.params_ = final
    self.history_ = [-h for h in history]
    return self

projected_kernel_matrix 

projected_kernel_matrix(
    feature_map: FeatureMap,
    X: NDArray[Any],
    gamma: float = 1.0,
    backend: BackendLike = None,
) -> NDArray[Any]

Projected quantum kernel — the standard answer to exponential concentration.

Instead of a global fidelity, compare the one-qubit reduced density matrices:

.. math:: k(x, x') = \exp\left(-\gamma \sum_i |\rho_i(x) - \rho_i(x')|_F^2\right)

Global overlaps concentrate as the register widens — every pair of inputs ends up looking equally similar, and the kernel stops carrying information. Local reduced states do not, so this stays informative where the fidelity kernel has already collapsed (Huang et al. 2021).

Source code in src/qmlkit/kernels/models.py 
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
def projected_kernel_matrix(
    feature_map: FeatureMap,
    X: npt.NDArray[Any],
    gamma: float = 1.0,
    backend: BackendLike = None,
) -> npt.NDArray[Any]:
    r"""Projected quantum kernel — the standard answer to exponential concentration.

    Instead of a *global* fidelity, compare the **one-qubit reduced density
    matrices**:

    .. math::  k(x, x') = \exp\left(-\gamma \sum_i \|\rho_i(x) - \rho_i(x')\|_F^2\right)

    Global overlaps concentrate as the register widens — every pair of inputs ends
    up looking equally similar, and the kernel stops carrying information. Local
    reduced states do not, so this stays informative where the fidelity kernel has
    already collapsed (Huang et al. 2021).
    """
    from qmlkit.core.execute import statevector

    rows = np.atleast_2d(np.asarray(X, dtype=float))
    n = feature_map.n_qubits
    # one reduced density matrix per qubit per sample
    rdms = []
    for r in rows:
        psi = statevector(feature_map.build(r), backend=backend)
        rdms.append([reduced_dm(psi, [q], n) for q in range(n)])

    m = len(rows)
    out = np.ones((m, m))
    for i in range(m):
        for j in range(i + 1, m):
            dist = sum(float(np.linalg.norm(rdms[i][q] - rdms[j][q], "fro") ** 2) for q in range(n))
            out[i, j] = out[j, i] = float(np.exp(-gamma * dist))
    return out

rkhs_model 

rkhs_model(
    alphas: Sequence[float],
    anchors: NDArray[Any],
    x: Sequence[float],
    kernel: Any,
) -> float

f(x) = sum_i alpha_i k(x_i, x) — a kernel model is a weighted similarity sum.

Source code in src/qmlkit/kernels/models.py 
306
307
308
309
310
311
312
313
314
315
def rkhs_model(
    alphas: Sequence[float],
    anchors: npt.NDArray[Any],
    x: Sequence[float],
    kernel: Any,
) -> float:
    """``f(x) = sum_i alpha_i k(x_i, x)`` — a kernel model is a weighted similarity sum."""
    a = np.asarray(alphas, dtype=float)
    pts = np.atleast_2d(np.asarray(anchors, dtype=float))
    return float(sum(ai * kernel(p, x) for ai, p in zip(a, pts, strict=True)))