# Quantum Reservoir Computing

Quantum reservoir computing (QRC) uses a quantum circuit as a fixed dynamical system, called the **reservoir**, that is never trained. Only a classical model on top of the features it produces, called the **readout**, is trained, and it is inexpensive to train. The weights of the reservoir are sampled once and then frozen, so all of the learning happens in the readout.

This page covers `eigenkit.reservoir`, the part of the library that uses only NumPy and needs no other framework. It contains the driver, the state policies, the feature extractors and a ridge-regression readout. Use these building blocks when you want full control, for example to process the features in your own way or to use your own readout from scikit-learn or PyTorch. The scikit-learn-style estimators built on them, `QuantumReservoirClassifier` and `QuantumReservoirRegressor`, are described on the [quantum machine learning (QML)](https://docs.kit.eigenorbit.com/user-guide/qml/index.md) page.

Everything relies on `Circuit.statevector(weights, inputs, initial_state=)`. The state is carried from one timestep to the next by reading the statevector out and injecting it back as the initial state of the following step.

## The driver

`run_reservoir(circuit, policy, sequences, weights)` evolves the circuit over a batch of sequences under a policy, and returns the trajectory of raw statevectors. The driver does not compute features or fit a readout, so the three layers can be combined freely: the driver first, then a feature extractor, then a readout.

```python
import numpy as np
from eigenkit import Circuit, Encoding, Variational, all_z
from eigenkit.reservoir import (
    run_reservoir,
    LeakyCarry,
    expectation_features,
    LinearReadout,
)

circuit = Circuit(
    blocks=[Encoding(in_features=3, rotation="RY"), Variational(depth=3, rotation="RY")],
    observables=all_z(3),
    provider="numpy",
)
weights = np.random.default_rng(0).standard_normal(circuit.spec.n_weights)
sequences = np.random.default_rng(1).uniform(-1, 1, size=(8, 40, 3))

states = run_reservoir(
    circuit, LeakyCarry(alpha=0.5), sequences, weights
)  # (8, 40, 8) complex
feats = expectation_features(states, "Z0 + Z1 + 0.5 X0 X1")  # (8, 40, 1)
readout = LinearReadout(ridge=1e-4)
```

## State policies

A *policy* decides, at each timestep, which state the circuit starts from and which inputs drive it. Two policies carry the quantum state from one timestep to the next, so the timesteps run one after another. The third policy carries no state, so all timesteps can run in parallel.

| Policy       | Constructor                     | What it does                                                            | Memory                                                        | Execution              |
| ------------ | ------------------------------- | ----------------------------------------------------------------------- | ------------------------------------------------------------- | ---------------------- |
| `Stateless`  | `Stateless(window=1)`           | Starts each step from \`                                                | 0…0⟩`and encodes the last`window\` inputs as rotation angles. | classical (the window) |
| `Carry`      | `Carry()`                       | Applies the circuit to the previous state, `ψ_t = U(x_t) ψ_{t-1}`.      | quantum state                                                 | sequential             |
| `LeakyCarry` | `LeakyCarry(alpha, reset=None)` | Mixes the carried state with a reset state and renormalises the result. | quantum, fading                                               | sequential             |

With the carry policies, the circuit takes `n_inputs = F` inputs per step, where `F` is the number of features per step. With `Stateless`, it takes the whole window, so the circuit must encode `n_inputs = window * F`. The driver checks this and raises a clear error otherwise.

Fading memory matters

A `Carry` policy preserves the norm of the state, but it tends to mix the history into a form that a linear readout cannot recover, so it often has little usable memory. A memory that fades over time is what makes a reservoir work in practice, so `LeakyCarry` with `alpha` between `0.3` and `0.7` is usually the right starting point for temporal tasks. The [memory-capacity example](https://docs.kit.eigenorbit.com/examples/reservoir-memory/index.md) shows this.

## Feature extractors

The readout trains on real-valued features extracted from the trajectory of statevectors. Each extractor takes the trajectory returned by the driver and returns real-valued features for each timestep.

| Extractor                                   | Features per step  | Notes                                                                                                                                                                                                         |
| ------------------------------------------- | ------------------ | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| `pauli_z_features(states)`                  | `n_qubits`         | The expectation value of `Z` on each qubit. This is the usual default in QRC; it is inexpensive and can be measured on hardware.                                                                              |
| `probability_features(states)`              | `2**n_qubits`      | The measurement distribution over all basis states. It carries more information, and it is quadratic in the amplitudes.                                                                                       |
| `amplitude_features(states)`                | `2 * 2**n_qubits`  | The real and imaginary parts of the state. Only a simulator can provide them; they cannot be measured on hardware.                                                                                            |
| `expectation_features(states, observables)` | one per observable | The expectation values of any Pauli observables. Pass the same `Observable` or `PauliExpr` objects, or Hamiltonian strings, that the rest of the library uses. A single Hamiltonian string is one observable. |

## Readout

`LinearReadout(ridge=1e-06)` performs ridge regression in closed form, with a bias term and no iterative optimiser. It is only a convenience: the features are a plain NumPy array, so any scikit-learn or PyTorch model can replace it.

## Backends and limitations

- Quantum reservoir computing runs on any backend that provides a statevector **and** accepts an injected initial state: NumPy, PennyLane, the Qiskit statevector simulator and the CUDA-Q simulators. Backends without access to the statevector, such as the Braket cloud services and Qiskit Runtime, raise a clear `NotImplementedError`. The default `provider="numpy"` needs no extra dependencies.
- The carry policies use state injection, which cannot be differentiated, so the reservoir is **frozen** and only the classical readout is trained. This is the standard setup in QRC, and it avoids computing gradients through the carried state.
- `LeakyCarry` approximates a dissipative reservoir with pure states. True dissipation needs density matrices, which the pure-state backends do not provide. The renormalised mixing used here approximates that behaviour, but it is not a physical quantum channel.

See the [short-term memory example](https://docs.kit.eigenorbit.com/examples/reservoir-memory/index.md) for a worked benchmark that compares the policies.

## API reference

### `eigenkit.reservoir`

Quantum reservoir computing (QRC).

A circuit with fixed weights, usually random, is driven over input sequences. The quantum state is carried from one timestep to the next, and classical features are read out from it.

The reservoir circuit is never trained; only a classical readout is. The driver, the policies, the feature extractors and the ridge-regression readout use only NumPy. Like `optimization` and `chemistry`, this package belongs to the core of the library, so it needs neither PyTorch nor a backend package.

Three policies cover the common variants of QRC:

- Stateless starts every step from `|0…0⟩` and encodes a sliding window of inputs. All timesteps run in parallel, in one batched evaluation.
- Carry carries the state from one step to the next with the circuit's unitary: `ψ_t = U(x_t) ψ_{t-1}`.
- LeakyCarry mixes the carried state with a reset state, and renormalises the result.

Example

```python
from eigenkit import Circuit
from eigenkit.reservoir import (
    run_reservoir, Carry, expectation_features, LinearReadout,
)

states = run_reservoir(circuit, Carry(), sequences, weights)  # (B, T, 2**n)
feats = expectation_features(states, "Z0 + Z1 + 0.5 X0 X1")   # (B, T, 1)
readout = LinearReadout(ridge=1e-4).fit(feats.reshape(B * T, -1), targets)
```

#### `Carry`

Policy that applies the circuit to the carried state: `ψ_t = U(x_t) ψ_{t-1}`.

The quantum state is the memory of the reservoir, and the first step starts from `|0…0⟩`. This policy is equivalent to `LeakyCarry(alpha=1.0)`.

##### `advance(x_t, evolved, state)`

Build the state for the next timestep from the result of this step.

Parameters:

| Name      | Type             | Description                                  | Default    |
| --------- | ---------------- | -------------------------------------------- | ---------- |
| `x_t`     | `ndarray`        | The input of this timestep, of shape (B, F). | *required* |
| `evolved` | `ndarray`        | The statevector that the circuit produced.   | *required* |
| `state`   | `ReservoirState` | The state that went into this step.          | *required* |

Returns:

| Type             | Description                      |
| ---------------- | -------------------------------- |
| `ReservoirState` | The state for the next timestep. |

##### `initial(batch)`

Return the state carried into the first timestep.

Parameters:

| Name    | Type  | Description                                   | Default    |
| ------- | ----- | --------------------------------------------- | ---------- |
| `batch` | `int` | The number of sequences that run in parallel. | *required* |

Returns:

| Type             | Description                 |
| ---------------- | --------------------------- |
| `ReservoirState` | The initial ReservoirState. |

##### `inject(x_t, state)`

Turn the input of this timestep into the arguments of the circuit call.

Parameters:

| Name    | Type             | Description                                    | Default    |
| ------- | ---------------- | ---------------------------------------------- | ---------- |
| `x_t`   | `ndarray`        | The input of this timestep, of shape (B, F).   | *required* |
| `state` | `ReservoirState` | The state carried over from the previous step. | *required* |

Returns:

| Name            | Type      | Description                                     |
| --------------- | --------- | ----------------------------------------------- |
| `initial_state` | \`ndarray | None\`                                          |
| `gate_inputs`   | `ndarray` | The values to pass to the gates of the circuit. |

#### `LeakyCarry`

Policy that mixes the carried state with a reset state and renormalises it.

It approximates a leaky, dissipative reservoir with pure states: `ψ_t = normalise(α U(x_t) ψ_{t-1} + (1 − α) reset)`.

`reset` defaults to `|0…0⟩` and should be a normalised statevector; this is not checked. With `alpha=1`, the policy is exactly Carry. This is not a true dissipative channel, which would need density matrices that the pure-state backends do not provide. It is an approximation that uses pure states only.

Attributes:

| Name    | Type      | Description                                                                                                                    |
| ------- | --------- | ------------------------------------------------------------------------------------------------------------------------------ |
| `alpha` | `float`   | The mixing weight of the carried state, in [0, 1]. Lower values forget faster, and 1 keeps everything, which is exactly Carry. |
| `reset` | \`ndarray | None\`                                                                                                                         |

##### `advance(x_t, evolved, state)`

Build the state for the next timestep from the result of this step.

Parameters:

| Name      | Type             | Description                                  | Default    |
| --------- | ---------------- | -------------------------------------------- | ---------- |
| `x_t`     | `ndarray`        | The input of this timestep, of shape (B, F). | *required* |
| `evolved` | `ndarray`        | The statevector that the circuit produced.   | *required* |
| `state`   | `ReservoirState` | The state that went into this step.          | *required* |

Returns:

| Type             | Description                      |
| ---------------- | -------------------------------- |
| `ReservoirState` | The state for the next timestep. |

##### `initial(batch)`

Return the state carried into the first timestep.

Parameters:

| Name    | Type  | Description                                   | Default    |
| ------- | ----- | --------------------------------------------- | ---------- |
| `batch` | `int` | The number of sequences that run in parallel. | *required* |

Returns:

| Type             | Description                 |
| ---------------- | --------------------------- |
| `ReservoirState` | The initial ReservoirState. |

##### `inject(x_t, state)`

Turn the input of this timestep into the arguments of the circuit call.

Parameters:

| Name    | Type             | Description                                    | Default    |
| ------- | ---------------- | ---------------------------------------------- | ---------- |
| `x_t`   | `ndarray`        | The input of this timestep, of shape (B, F).   | *required* |
| `state` | `ReservoirState` | The state carried over from the previous step. | *required* |

Returns:

| Name            | Type      | Description                                     |
| --------------- | --------- | ----------------------------------------------- |
| `initial_state` | \`ndarray | None\`                                          |
| `gate_inputs`   | `ndarray` | The values to pass to the gates of the circuit. |

#### `LinearReadout`

Ridge regression in closed form, with a bias term.

It solves `W = (XᵀX + λI)⁻¹ XᵀY`, and the bias column is not regularised. `Y` can have one output or several. `ridge` is the strength `λ` of the L2 penalty, which adds `λ‖W‖²` to the squared error of the fit. A small positive value also keeps the solution numerically stable when the reservoir features are collinear.

Parameters:

| Name    | Type    | Description                                                                                                                                                                   | Default |
| ------- | ------- | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ------- |
| `ridge` | `float` | The strength λ of the L2 penalty. Raise it when the readout overfits or the feature matrix is nearly singular, and lower it when the readout underfits. The default is 1e-06. | `1e-06` |

##### `fit(X, y)`

Fit the readout to `X` and `y`.

Parameters:

| Name | Type      | Description                                                                               | Default    |
| ---- | --------- | ----------------------------------------------------------------------------------------- | ---------- |
| `X`  | `ndarray` | The samples, of shape (n_samples, n_features).                                            | *required* |
| `y`  | `ndarray` | The targets, of shape (n_samples,) for one output, or (n_samples, n_outputs) for several. | *required* |

Returns:

| Type              | Description                                 |
| ----------------- | ------------------------------------------- |
| `'LinearReadout'` | This readout, so that calls can be chained. |

##### `predict(X)`

Predict the target for every sample in `X`.

Parameters:

| Name | Type      | Description                                    | Default    |
| ---- | --------- | ---------------------------------------------- | ---------- |
| `X`  | `ndarray` | The samples, of shape (n_samples, n_features). | *required* |

Returns:

| Type      | Description                                                                                           |
| --------- | ----------------------------------------------------------------------------------------------------- |
| `ndarray` | One prediction per sample, with one column per output when the readout was fitted on several outputs. |

#### `ReservoirPolicy`

Bases: `Protocol`

A reservoir policy that carries state from one timestep to the next.

Parallel policies, which carry no state, do not implement this protocol. They provide `expand_inputs(sequences)`, which returns the inputs of every timestep with shape `(B, T, n_inputs)`, and the driver evaluates all timesteps at once. See Stateless.

##### `advance(x_t, evolved, state)`

Build the state for the next timestep from the result of this step.

Parameters:

| Name      | Type             | Description                                  | Default    |
| --------- | ---------------- | -------------------------------------------- | ---------- |
| `x_t`     | `ndarray`        | The input of this timestep, of shape (B, F). | *required* |
| `evolved` | `ndarray`        | The statevector that the circuit produced.   | *required* |
| `state`   | `ReservoirState` | The state that went into this step.          | *required* |

Returns:

| Type             | Description                      |
| ---------------- | -------------------------------- |
| `ReservoirState` | The state for the next timestep. |

##### `initial(batch)`

Return the state carried into the first timestep.

Parameters:

| Name    | Type  | Description                                   | Default    |
| ------- | ----- | --------------------------------------------- | ---------- |
| `batch` | `int` | The number of sequences that run in parallel. | *required* |

Returns:

| Type             | Description                 |
| ---------------- | --------------------------- |
| `ReservoirState` | The initial ReservoirState. |

##### `inject(x_t, state)`

Turn the input of this timestep into the arguments of the circuit call.

Parameters:

| Name    | Type             | Description                                    | Default    |
| ------- | ---------------- | ---------------------------------------------- | ---------- |
| `x_t`   | `ndarray`        | The input of this timestep, of shape (B, F).   | *required* |
| `state` | `ReservoirState` | The state carried over from the previous step. | *required* |

Returns:

| Name            | Type      | Description                                     |
| --------------- | --------- | ----------------------------------------------- |
| `initial_state` | \`ndarray | None\`                                          |
| `gate_inputs`   | `ndarray` | The values to pass to the gates of the circuit. |

#### `ReservoirState`

The state that sequential policies carry from one timestep to the next.

`vector` is the `(B, 2**n)` statevector to inject at the next step, and `None` means that the step starts from `|0…0⟩`. `history` holds the recent raw inputs, of shape `(B, k, F)`, for a sequential policy that encodes a window of inputs.

Attributes:

| Name      | Type      | Description |
| --------- | --------- | ----------- |
| `vector`  | \`ndarray | None\`      |
| `history` | \`ndarray | None\`      |

#### `Stateless`

Parallel policy that carries no state.

Every step starts from `|0…0⟩`, and the last `window` inputs are encoded as rotation angles of the circuit's gates. The memory is purely classical (the window), and the quantum state carries nothing from one step to the next. Every timestep is therefore independent, and the driver evaluates all of them in one batched run. With `window=1`, the policy has no memory: it encodes only the current input.

The circuit's `n_inputs` must equal `window * F`, where `F` is the number of features per step, because each step passes the whole flattened window to the encoding gates.

Attributes:

| Name     | Type  | Description                                                                                                                                                                                        |
| -------- | ----- | -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| `window` | `int` | The number of consecutive inputs to encode at each step. 1 means no memory: only the current input is encoded. Whatever the value, the circuit's n_inputs must equal window * F. The default is 1. |

##### `expand_inputs(sequences)`

Return the gate inputs of every timestep, with shape `(B, T, window * F)`.

Step `t` encodes the inputs `t - window + 1` to `t`. The first `window - 1` steps are padded with zeros at the front. No quantum state is involved, so this is only a reshaping of the input sequence.

Parameters:

| Name        | Type      | Description                                | Default    |
| ----------- | --------- | ------------------------------------------ | ---------- |
| `sequences` | `ndarray` | The whole input batch, of shape (B, T, F). | *required* |

Returns:

| Type      | Description                                                      |
| --------- | ---------------------------------------------------------------- |
| `ndarray` | The gate inputs for every timestep, of shape (B, T, window * F). |

#### `amplitude_features(states)`

The raw statevector as real features: its real and imaginary parts, concatenated.

`states` is a complex array of shape `(B, T, 2**n)`, and the result is a real array of shape `(B, T, 2 * 2**n)`. It exposes the full state, which only a simulator can provide; it cannot be measured on hardware.

Parameters:

| Name     | Type      | Description                                              | Default    |
| -------- | --------- | -------------------------------------------------------- | ---------- |
| `states` | `ndarray` | The trajectory of statevectors, of shape (B, T, 2\*\*n). | *required* |

Returns:

| Type      | Description                                                    |
| --------- | -------------------------------------------------------------- |
| `ndarray` | The real and imaginary parts, concatenated, for each timestep. |

#### `expectation_features(states, observables)`

Expectation values `⟨ψ|H_k|ψ⟩` of each observable, over the whole trajectory.

`states` is a complex array of shape `(B, T, 2**n)`, and the result is a real array of shape `(B, T, n_obs)`. `observables` can be an `Observable`, a `PauliExpr` or a Hamiltonian string, or a sequence of these. These are the same types that the rest of the library uses, so measuring an arbitrary Hamiltonian, as QRC experiments often require, needs no manual handling of Pauli terms. A single Hamiltonian string is one observable, and gives one feature.

Parameters:

| Name          | Type                   | Description                                              | Default    |
| ------------- | ---------------------- | -------------------------------------------------------- | ---------- |
| `states`      | `ndarray`              | The trajectory of statevectors, of shape (B, T, 2\*\*n). | *required* |
| `observables` | `Sequence[Observable]` | The observables to measure.                              | *required* |

Returns:

| Type      | Description                                  |
| --------- | -------------------------------------------- |
| `ndarray` | One feature per observable and per timestep. |

#### `pauli_z_features(states)`

The expectation value `⟨Z_q⟩` on each qubit, the usual default readout in QRC.

`states` has shape `(B, T, 2**n)`, and the result has shape `(B, T, n)`. The values are computed with expectation_features(), so both functions share the same tested measurement code.

Parameters:

| Name     | Type      | Description                                              | Default    |
| -------- | --------- | -------------------------------------------------------- | ---------- |
| `states` | `ndarray` | The trajectory of statevectors, of shape (B, T, 2\*\*n). | *required* |

Returns:

| Type      | Description                             |
| --------- | --------------------------------------- |
| `ndarray` | One feature per qubit and per timestep. |

#### `probability_features(states)`

The measurement distribution over all `2**n` basis states, `|ψ|²`.

`states` is a complex array of shape `(B, T, 2**n)`, and the result is a real array of the same shape.

Parameters:

| Name     | Type      | Description                                              | Default    |
| -------- | --------- | -------------------------------------------------------- | ---------- |
| `states` | `ndarray` | The trajectory of statevectors, of shape (B, T, 2\*\*n). | *required* |

Returns:

| Type      | Description                                                |
| --------- | ---------------------------------------------------------- |
| `ndarray` | The distribution over the basis states, for each timestep. |

#### `run_reservoir(circuit, policy, sequences, weights)`

Evolve `circuit` over `sequences` under `policy`, with frozen `weights`.

The result is the trajectory of statevectors, a complex array of shape `(B, T, 2**n)`.

`sequences` has shape `(B, T, F)`: `B` independent series, `T` timesteps and `F` features per step. `weights` is the fixed 1-D parameter vector of the reservoir, shared by the whole batch. Features are extracted separately from the returned trajectory (see the feature extractors), so the same run can produce expectation, probability or amplitude features without evolving the circuit again.

Parameters:

| Name        | Type      | Description                                         | Default    |
| ----------- | --------- | --------------------------------------------------- | ---------- |
| `circuit`   | `Circuit` | The reservoir circuit, run with frozen weights.     | *required* |
| `policy`    | `Policy`  | How the state passes from one timestep to the next. | *required* |
| `sequences` | `ndarray` | The inputs, of shape (B, T, F).                     | *required* |
| `weights`   | `ndarray` | The frozen weights of the reservoir.                | *required* |

Returns:

| Type      | Description                                              |
| --------- | -------------------------------------------------------- |
| `ndarray` | The trajectory of statevectors, of shape (B, T, 2\*\*n). |