Nonlinear regression with a quantum extreme learning machine

Quantum Machine Learning, worked example.

Open In Colab

Nonlinear regression with a quantum extreme learning machine

A fixed quantum system turns a single number \(x\) into a bank of nonlinear features. Only a classical linear readout is trained. This notebook builds such a quantum extreme learning machine (QELM) from scratch with exact statevectors (NumPy and SciPy only, no quantum-computing library) and uses it to regress a multi-frequency target function.

Motivation: the extreme learning machine idea

A classical extreme learning machine (ELM) is a neural network with one hidden layer whose input-to-hidden weights are set randomly and never trained. The hidden layer applies a fixed nonlinear map \(\phi\) to the input, producing a high-dimensional feature vector, and only the hidden-to-output weights are fit, by linear least squares. The bet is that a rich enough random feature map already contains, in linear combination, the function one wants; learning then reduces to a single convex problem with no backpropagation and no local minima.

The quantum version replaces the random hidden layer by a fixed quantum process. A classical input \(x\) is encoded into a quantum state, that state is processed by a fixed, untrained unitary, and a set of observables is measured. The vector of expectation values is the feature vector \(\boldsymbol{f}(x)\). As in the classical ELM, only a linear readout \(\boldsymbol{w}\) is trained:

\[ y(x) \;=\; \boldsymbol{w}^{\mathsf T}\,\boldsymbol{f}(x), \qquad f_\alpha(x) \;=\; \langle \hat{O}_\alpha \rangle_{\rho_{\mathrm{out}}(x)} . \]

Two features distinguish the QELM from a quantum reservoir computer (QRC). The QELM is static: each input \(x\) is processed independently, with no memory of previous inputs. And its expressive power comes from the encoding and the fixed scrambling dynamics, not from any temporal recurrence.

Problem statement

Regress the multi-frequency target

\[ f(x) \;=\; \sin(3x) \;+\; \tfrac{1}{2}\sin(7x), \qquad x \in [0, 2\pi], \]

from a handful of training samples, using a QELM whose reservoir is a fixed, nonintegrable spin chain.

Protocol for one input \(x\).

  1. Start from \(|0\rangle^{\otimes N}\) on a small register (\(N=5\)).
  2. Repeat \(L\) times: apply a fixed scrambling unitary \(\hat U = e^{-i\hat H t}\) with \(\hat H\) the mixed-field Ising Hamiltonian, then a data-encoding layer \(\hat S(x) = \prod_i e^{-i x \hat Z_i / 2}\). This interleaving is data re-uploading.
  3. Apply one final scramble, so the encoded phases become measurable coherences.
  4. Read a fixed bank of low-weight Pauli expectation values as \(\boldsymbol{f}(x)\).

Only the ridge-regression readout \(\boldsymbol{w}\) is trained. We report the held-out \(R^2\) and mean squared error, and we show that the number of re-uploads \(L\) controls which Fourier frequencies the reservoir can represent.

Code 
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import expm

np.set_printoptions(precision=4, suppress=True)

# Single-qubit Pauli matrices.
I2 = np.eye(2, dtype=complex)
X = np.array([[0, 1], [1, 0]], dtype=complex)
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)


def kron_op(single, site, N):
    "Embed a single-qubit operator at `site` into the N-qubit space."
    op = np.array([[1.0]], dtype=complex)
    for j in range(N):
        op = np.kron(op, single if j == site else I2)
    return op

The reservoir: a fixed scrambling Hamiltonian

The reservoir is the nonintegrable mixed-field Ising chain

\[ \hat H \;=\; J\sum_i \hat Z_i \hat Z_{i+1} \;+\; h_x \sum_i \hat X_i \;+\; h_z \sum_i \hat Z_i . \]

The transverse field \(h_x\) alone would leave the model Jordan-Wigner integrable; the longitudinal field \(h_z\) breaks integrability, so the dynamics genuinely scramble and mix low-weight operators into high-weight ones. This is what supplies a rich, nonlinear feature basis. The unitary \(\hat U = e^{-i\hat H t}\) is computed once and never trained.

Code 
def mixed_field_ising(N, J, hx, hz):
    "Nonintegrable mixed-field Ising Hamiltonian on N qubits."
    H = np.zeros((2 ** N, 2 ** N), dtype=complex)
    for i in range(N - 1):
        H += J * (kron_op(Z, i, N) @ kron_op(Z, i + 1, N))
    for i in range(N):
        H += hx * kron_op(X, i, N) + hz * kron_op(Z, i, N)
    return H


N = 5                         # small register: exact statevectors, dim 2^5 = 32
J, hx, hz = 1.0, 0.9, 0.8     # nonintegrable point
t = 0.6                       # per-layer scrambling time

H = mixed_field_ising(N, J, hx, hz)
U = expm(-1j * H * t)         # fixed reservoir unitary, computed once
print("Hilbert space dimension:", 2 ** N)
print("U is unitary:", np.allclose(U @ U.conj().T, np.eye(2 ** N)))
Hilbert space dimension: 32
U is unitary: True

Encoding and the Fourier picture

The input is written in by the diagonal encoding layer \(\hat S(x) = \prod_i e^{-i x \hat Z_i / 2}\). Acting on a statevector it multiplies each computational-basis amplitude \(b\) by a phase \(e^{-i x g_b}\), where \(g_b = \tfrac12\sum_i (1 - 2\,b_i)\) is the total \(\hat Z/2\) eigenvalue of that basis state.

Every measured feature is therefore a finite Fourier series in \(x\). The accessible frequencies are the differences of the encoding generator’s eigenvalues, and with \(L\) re-uploads these differences add up: the reachable frequencies extend to \(L\cdot N\). This is the Schuld-Sweke-Meyer picture. With \(N=5\) a single upload reaches frequency \(5\), which cannot represent the \(\sin(7x)\) component of the target; two uploads reach frequency \(10\) and can. Re-uploads, not extra training, buy the expressivity.

Code 
def encoding_phases(N):
    "Per-basis-state phase generator g_b for S(x) = prod_i exp(-i x Z_i / 2)."
    b = np.arange(1 << N)
    pop = np.array([bin(i).count("1") for i in range(1 << N)])
    return (N - 2 * pop) / 2.0     # g_b: eigenvalue differences are integers 0..N


g = encoding_phases(N)
print("encoding generator eigenvalues g_b:", np.unique(g))
print("max reachable frequency at L uploads = L * N =  L *", N)
encoding generator eigenvalues g_b: [-2.5 -1.5 -0.5  0.5  1.5  2.5]
max reachable frequency at L uploads = L * N =  L * 5

Measurement: the readout observables

The feature vector collects low-weight Pauli expectation values: all single-site \(\{\hat X_i, \hat Y_i, \hat Z_i\}\), the nearest-neighbour \(\{\hat X_i\hat X_{i+1}, \hat Y_i\hat Y_{i+1}, \hat Z_i\hat Z_{i+1}\}\), and the all-pairs \(\{\hat Z_i \hat Z_j\}\). For a pure state \(\psi\) any Pauli string \(\hat P = i^p \hat X^{x}\hat Z^{z}\) has

\[ \langle \hat P\rangle = \mathrm{Re}\sum_b \overline{\psi_b}\, W_b\, \psi_{b\oplus x}, \]

an \(\mathcal{O}(2^N)\) contraction that avoids building the full matrix. We verify it against the brute-force expectation below.

Code 
_POP = np.array([bin(i).count("1") for i in range(1 << 16)])
_PH = np.array([1, 1j, -1, -1j])   # i^p


def pauli_readout(N):
    "Column-index and weight arrays for the low-weight Pauli readout bank."
    bit = lambda i: 1 << (N - 1 - i)
    specs, names = [], []
    for i in range(N):
        specs += [(bit(i), 0, 0), (bit(i), bit(i), 1), (0, bit(i), 0)]
        names += [f"X{i}", f"Y{i}", f"Z{i}"]
    for i in range(N - 1):
        m = bit(i) | bit(i + 1)
        specs += [(m, 0, 0), (m, m, 2), (0, m, 0)]
        names += [f"X{i}X{i+1}", f"Y{i}Y{i+1}", f"Z{i}Z{i+1}"]
    for i in range(N):
        for j in range(i + 2, N):
            m = bit(i) | bit(j)
            specs += [(0, m, 0)]
            names += [f"Z{i}Z{j}"]
    b = np.arange(1 << N)
    COLS = np.array([b ^ x for (x, z, p) in specs])
    W = np.array([_PH[p] * (1 - 2 * (_POP[z & b] & 1)) for (x, z, p) in specs],
                 dtype=complex)
    return COLS, W, names


COLS, W, names = pauli_readout(N)
print(f"{len(names)} readout observables:", names)

# Sanity check: the fast contraction reproduces the brute-force <psi|P|psi>.
rng = np.random.default_rng(0)
psi = rng.normal(size=1 << N) + 1j * rng.normal(size=1 << N)
psi /= np.linalg.norm(psi)
fast = np.real(np.conj(psi) * W * psi[COLS]).sum(axis=1)
X0 = kron_op(X, 0, N)
Z0Z1 = kron_op(Z, 0, N) @ kron_op(Z, 1, N)
print("check <X0>:", np.allclose(fast[names.index("X0")], (psi.conj() @ X0 @ psi).real))
print("check <Z0Z1>:", np.allclose(fast[names.index("Z0Z1")], (psi.conj() @ Z0Z1 @ psi).real))
33 readout observables: ['X0', 'Y0', 'Z0', 'X1', 'Y1', 'Z1', 'X2', 'Y2', 'Z2', 'X3', 'Y3', 'Z3', 'X4', 'Y4', 'Z4', 'X0X1', 'Y0Y1', 'Z0Z1', 'X1X2', 'Y1Y2', 'Z1Z2', 'X2X3', 'Y2Y3', 'Z2Z3', 'X3X4', 'Y3Y4', 'Z3Z4', 'Z0Z2', 'Z0Z3', 'Z0Z4', 'Z1Z3', 'Z1Z4', 'Z2Z4']
check <X0>: True
check <Z0Z1>: True

The QELM feature map

Now assemble the single-shot map \(x \mapsto \boldsymbol{f}(x)\). Nothing here is trained: the reservoir unitary \(\hat U\), the encoding, and the readout are all fixed. The only nonlinearity the readout ever sees is manufactured by this quantum pipeline.

Code 
def qelm_features(xs, N, U, g, COLS, W, L):
    "Feature matrix F (len(xs) x d): the fixed encode->scramble QELM at each x."
    F = np.zeros((len(xs), COLS.shape[0]))
    for n, x in enumerate(xs):
        psi = np.zeros(1 << N, dtype=complex)
        psi[0] = 1.0
        for _ in range(L):
            psi = U @ psi              # fixed scrambling
            psi *= np.exp(-1j * x * g) # data-encoding layer S(x)
        psi = U @ psi                  # final scramble -> measurable coherences
        F[n] = np.real(np.conj(psi) * W * psi[COLS]).sum(axis=1)
    return F


def target(x):
    return np.sin(3 * x) + 0.5 * np.sin(7 * x)

Visualising the feature bank

A single scalar \(x\) is turned into a family of nonlinear, multi-frequency curves. These are the columns of the feature matrix, plotted against \(x\). The linear readout will fit the target as a weighted sum of exactly these curves.

Code 
x_dense = np.linspace(0, 2 * np.pi, 400)
F_dense = qelm_features(x_dense, N, U, g, COLS, W, L=2)

plt.figure(figsize=(7, 4))
for nm in ["X0", "Y2", "Z4", "Z0Z1", "Y1"]:
    plt.plot(x_dense, F_dense[:, names.index(nm)], lw=1.4, label=fr"$\langle {nm} \rangle$")
plt.xlabel("input $x$"); plt.ylabel(r"feature $\langle O_\alpha\rangle(x)$")
plt.xticks([0, np.pi, 2 * np.pi], ["0", r"$\pi$", r"$2\pi$"])
plt.legend(ncol=2, fontsize=9); plt.title("Untrained reservoir features")
plt.tight_layout(); plt.show()

Training the linear readout

We standardise the features on the training split, append a bias term, and solve the ridge-regression normal equations. This is the only fit in the whole pipeline, a single convex linear-algebra step with no variational loop.

Code 
def ridge_fit(F_tr, y_tr, lam):
    mu, sd = F_tr.mean(0), F_tr.std(0) + 1e-12
    Z = np.hstack([(F_tr - mu) / sd, np.ones((len(F_tr), 1))])
    w = np.linalg.solve(Z.T @ Z + lam * np.eye(Z.shape[1]), Z.T @ y_tr)
    return w, mu, sd


def ridge_predict(F, w, mu, sd):
    Z = np.hstack([(F - mu) / sd, np.ones((len(F), 1))])
    return Z @ w


def r2_score(y, yhat):
    return 1.0 - np.sum((y - yhat) ** 2) / np.sum((y - y.mean()) ** 2)


L_main, lam = 2, 1e-8
x_train = np.linspace(0, 2 * np.pi, 80, endpoint=False)
x_test = np.linspace(0, 2 * np.pi, 400)
y_train, y_test = target(x_train), target(x_test)

F_tr = qelm_features(x_train, N, U, g, COLS, W, L_main)
F_te = qelm_features(x_test, N, U, g, COLS, W, L_main)
w, mu, sd = ridge_fit(F_tr, y_train, lam)
pred = ridge_predict(F_te, w, mu, sd)

r2 = r2_score(y_test, pred)
mse = np.mean((y_test - pred) ** 2)
print(f"held-out R^2 = {r2:.4f}")
print(f"held-out MSE = {mse:.3e}")
held-out R^2 = 1.0000
held-out MSE = 1.773e-06
Code 
plt.figure(figsize=(7, 4))
plt.plot(x_test, y_test, "k-", lw=2, label="target $f(x)$")
plt.plot(x_test, pred, "r--", lw=1.6, label="QELM readout")
plt.scatter(x_train, y_train, s=18, color="#1f6fc4", edgecolors="white",
            linewidths=0.5, zorder=3, label="training samples")
plt.xlabel("input $x$"); plt.ylabel("$f(x)$")
plt.xticks([0, np.pi, 2 * np.pi], ["0", r"$\pi$", r"$2\pi$"])
plt.text(0.03, 0.05, f"test $R^2 = {r2:.3f}$", transform=plt.gca().transAxes,
         bbox=dict(boxstyle="round", fc="white", ec="0.7"))
plt.legend(loc="upper right", fontsize=9); plt.title("QELM regression")
plt.tight_layout(); plt.show()

Expressivity grows with re-uploads

The pedagogical payoff. Sweep the number of re-uploads \(L\) at a weak and a strong per-layer scrambling time. With one upload the reachable frequencies stop at \(N=5\), so the readout can only fit the \(\sin(3x)\) part and \(R^2\) is capped near \(0.80\), exactly the variance fraction of that component. Adding a second upload lifts the frequency ceiling past \(7\) and the fit becomes essentially exact, provided the scrambling is strong enough; a weakly scrambling reservoir needs an extra upload to build up the same high-frequency amplitude. Pushing \(L\) higher eventually over-scrambles and the feature signal degrades.

Code 
L_sweep = [1, 2, 3, 4]
for tev, label in [(0.15, "weak  "), (0.6, "strong")]:
    Ut = expm(-1j * H * tev)
    r2s = []
    for L in L_sweep:
        Ftr = qelm_features(x_train, N, Ut, g, COLS, W, L)
        Fte = qelm_features(x_test, N, Ut, g, COLS, W, L)
        wl, mul, sdl = ridge_fit(Ftr, y_train, lam)
        r2s.append(r2_score(y_test, ridge_predict(Fte, wl, mul, sdl)))
    plt.plot(L_sweep, r2s, "-o", label=f"{label} (t={tev})")
    print(f"{label}  R^2(L=1..4) = {np.round(r2s, 3).tolist()}")

plt.axhline(1.0, color="0.6", ls=":", lw=1)
plt.axvspan(0.5, 1.5, color="0.85", alpha=0.5)
plt.xlabel("number of re-uploads $L$"); plt.ylabel("test $R^2$")
plt.xticks(L_sweep); plt.ylim(-0.05, 1.08)
plt.legend(title="per-layer scrambling", fontsize=9)
plt.title("Expressivity vs re-uploads"); plt.tight_layout(); plt.show()
weak    R^2(L=1..4) = [0.799, 0.8, 0.995, 0.964]
strong  R^2(L=1..4) = [0.8, 1.0, 0.993, 0.833]

Summary

  • A fixed encode-then-scramble quantum pipeline maps a scalar input \(x\) into a bank of nonlinear, multi-frequency features in a single shot. No temporal memory is involved: every input is processed independently.
  • Only a ridge-regression readout is trained, a single convex step. On the multi-frequency target the held-out fit is essentially exact (\(R^2 \approx 1\)).
  • Expressivity is set by the encoding spectrum and the number of re-uploads, not by training: with \(N\) qubits and \(L\) uploads the reachable Fourier frequencies extend to \(L\cdot N\). One upload cannot represent \(\sin(7x)\); two can. This is the data-reuploading mechanism, with scrambling supplying the cross-frequency amplitudes.
  • The scrambling strength is the second knob: stronger per-layer mixing reaches a given accuracy in fewer uploads, while too much mixing over-scrambles the features.