Photonic Circuits (Perceval)

Experimental

The photonics package is experimental. Its public interface may change, and it currently evaluates one sample and one scalar observable at a time. The roadmap lists what comes next.

eigenkit models linear-optics (photonic) circuits in a separate package, with its own model of computation. They are not a backend under the qubit Provider, because the two models of computation are different:

Qubit circuits (`core`)	Photonic circuits (`photonics`)
A statevector of size `2ⁿ`.	A Fock space, made of the configurations of photon numbers over the modes.
Gates act on qubits.	Beam splitters and phase shifters act on modes.
The readout is an expectation value of Pauli operators.	The readout is a distribution of photon numbers, or detection statistics.
The parameter-shift rule shifts by `±π/2`.	A photonic parameter-shift rule uses more terms.

Placing photonics under the qubit Provider would work only through dual-rail encoding. That encoding is probabilistic, because its two-qubit gates rely on post-selection, and it gives up everything that makes photonics interesting. So photonics has its own intermediate representation (IR) and its own executor. It reuses the parts of the core that do not depend on the model of computation: the Parameter expression tree, the separation between compiling and running, and the bridge to PyTorch through autograd.Function.

Install

pip install eigenkit[perceval]

This installs Perceval, the photonic software development kit of Quandela, under the MIT licence. It is not part of the all extra, because the Perceval simulator pins protobuf to an exact version. Install it explicitly, for example with eigenkit[all,perceval]. Importing eigenkit, or even eigenkit.photonics, never imports Perceval. Perceval is loaded only when you create or run a Perceval executor.

Quickstart: the Hong-Ou-Mandel dip

A circuit is a PhotonicCircuitSpec: a number of modes n_modes, an ordered tuple of components, and an input Fock state, which gives the number of photons in each mode. The angles of the components are the same Parameter expressions as in the qubit intermediate representation, or plain numbers.

import math
import numpy as np

from eigenkit.photonics import BeamSplitter, PercevalBackend, PhotonicCircuitSpec

# Two indistinguishable photons, one in each mode, into a balanced beam-splitter.
spec = PhotonicCircuitSpec(
    n_modes=2,
    n_weights=0,
    n_inputs=0,
    components=(BeamSplitter(modes=(0, 1), theta=math.pi / 2),),  # theta=pi/2 => 50:50
    input_state=(1, 1),
)

compiled = PercevalBackend().compile(spec)  # compile once
dist = compiled.run(np.array([]), np.array([]))  # {fock_tuple: probability}
print(dist)  # {(0, 2): 0.5, (2, 0): 0.5}: the photons bunch, and P(1, 1) = 0

The photons bunch: (1, 1) has probability 0. This result is not classical, because it comes from the permanent of the transfer matrix of the circuit. It is a good check that the executor computes real boson statistics.

The components available today are BeamSplitter(modes, theta, phi=0), PhaseShifter(mode, phi) and Permutation(perm). An angle can be a plain number (stored as a Constant), a WeightRef (trainable), an InputRef (data), or any combination with Scale, Sum and Multiply. This is exactly the Parameter tree of the qubit circuits, and it accepts the same values as the qubit CircuitBuilder.

Observables

The readouts are observables that are diagonal in the Fock basis, which are the ones that photonic quantum machine learning uses in practice. A FockProbability target must hold as many photons as the input state, because linear optics conserves the number of photons.

from eigenkit.photonics import FockProbability, PhotonNumber

w, x = np.array([]), np.array([])  # the circuit above has no weights and no inputs
compiled.expectation(w, x, PhotonNumber(mode=0))  # <n̂_0>, expected photons in mode 0
compiled.expectation(w, x, FockProbability(target=(2, 0)))  # P(detect exactly |2,0>)

Gradients: the photonic parameter-shift rule

The qubit rule with shifts of ±π/2 does not apply to linear optics. The generator of a phase shifter is the number operator, whose spectrum is made of integers. So an expectation that is diagonal in the Fock basis is a trigonometric polynomial in the phase, of degree at most N, where N is the total number of photons. Its exact derivative needs the general equidistant shift rule with 2N terms. The rule with ±π/2 is the case N = 1.

The angle theta of a beam splitter is generated by an operator of the special unitary group SU(2) on its two modes. The eigenvalues of this operator are integers or half-integers, depending on how many photons the two modes hold. In a circuit with two modes, this count is always N, and the same rule with 2N terms applies. With more modes, the expectation also has half-integer frequencies in theta, so the rule is applied in theta / 2 with order 2N, which takes 4N evaluations per beam-splitter angle. Both rules are exact. They reuse chain_factor from the qubit intermediate representation to map the derivatives with respect to the angles onto the weights and the inputs.

from eigenkit.core.ir import WeightRef
from eigenkit.photonics import PhaseShifter

# One photon through a Mach-Zehnder interferometer with a trainable phase.
spec = PhotonicCircuitSpec(
    n_modes=2,
    n_weights=1,
    n_inputs=0,
    components=(
        BeamSplitter(modes=(0, 1), theta=math.pi / 2),
        PhaseShifter(mode=0, phi=WeightRef(0)),
        BeamSplitter(modes=(0, 1), theta=math.pi / 2),
    ),
    input_state=(1, 0),
)
compiled = PercevalBackend().compile(spec)
w = np.array([0.3])
weight_jac, input_jac = compiled.jacobian(w, x, PhotonNumber(0), method="parameter_shift")

method="finite_difference" gives a gradient by central differences, as a check. With two or more photons, the naive rule with ±π/2 is wrong, and only the rule with more terms matches the finite differences (see tests/photonics/test_gradients.py). Pass need_input_grad=False, or need_weight_grad=False, to skip the shift evaluations for a side that you do not need. That side is then returned as None. PhotonicCircuit does this automatically for the inputs that do not require a gradient.

Training with PyTorch

PhotonicCircuit is an nn.Module that holds the weights as an nn.Parameter and passes the photonic gradient to PyTorch autograd. So loss.backward() and any torch.optim optimiser work as usual:

import torch
from eigenkit.photonics import FockProbability
from eigenkit.photonics import PhotonicCircuit

circuit = PhotonicCircuit(spec, FockProbability(target=(0, 1)), init_weights=[1.5])
opt = torch.optim.Adam(circuit.parameters(), lr=0.1)
for _ in range(120):
    opt.zero_grad()
    loss = (1.0 - circuit()) ** 2  # for example, send the photon to mode 1
    loss.backward()
    opt.step()

A complete worked classifier is in examples/photonics_classification.py. It encodes 2-D features as phases, trains a shallow interferometer with Adam, and reaches an accuracy of about 0.97.

Place the encodings of a photonic feature map so that they do not commute

When two features are encoded as phases on adjacent single modes of a one-photon state, they only change their relative phase. The model then sees x0 − x1 and cannot see x0 + x1. Separate the two encodings with a beam splitter, which does not commute with them, so that both features enter independently. This is the photonic version of ordering the rotations in a qubit feature map.

Reservoir computing

Photonics also supports reservoir computing, in the form of a quantum extreme learning machine. A fixed random interferometer is the reservoir. Each input is encoded in its phases and read out as photon-number statistics, which gives a high-dimensional nonlinear map of the features. Only a classical linear head is trained. No gradient through the optics is needed, only the forward path of the executor. It reuses the Stateless windowing policy of the qubit reservoir package without changes, and the readout can be the closed-form LinearReadout of the same package.

import numpy as np
from eigenkit.photonics.reservoir import random_reservoir_spec, run_photonic_reservoir
from eigenkit.reservoir import LinearReadout, Stateless

rng = np.random.default_rng(0)
sequences = rng.uniform(0, np.pi, size=(8, 10, 1))  # (B, T, F), one feature per step
targets = sequences[..., 0] * np.roll(sequences[..., 0], 1, axis=1)  # x_t * x_{t-1}

spec, weights = random_reservoir_spec(n_modes=4, n_inputs=2, input_state=(1, 1, 0, 0))
feats = run_photonic_reservoir(
    spec, weights, sequences, policy=Stateless(window=2)
)  # (B, T, n_fock)
head = LinearReadout(ridge=1e-6).fit(
    feats.reshape(-1, feats.shape[-1]), targets.reshape(-1)
)

A worked example is in examples/photonics_reservoir.py. It learns the bilinear memory task y_t = x_t·x_{t-1}, which a linear readout on the raw inputs cannot learn.

No memory yet

Only reservoirs without a state, or with a classical window, are supported. A photonic carry-over of the state, which would be a quantum memory across time steps, needs Fock states to be injected into the executor, and it is not supported yet.

Running on Quandela Cloud

Besides local simulation, a circuit can run on Quandela Cloud, on a real photonic quantum processing unit (QPU) or on a cloud simulator, through QuandelaCloudBackend. The results are sampled (counts of shots), so run returns an estimated Fock distribution, and gradients are not available, as with the qubit QPU backends. The backend follows the QPU conventions of the library. It issues a CostWarning on the first submission, submit() returns a handle with .result() and .status() without waiting, and dry_run_info() works offline.

import numpy as np
from eigenkit.photonics import QuandelaCloudBackend

backend = QuandelaCloudBackend("qpu:ascella", shots=20000)  # token from $QUANDELA_TOKEN
compiled = backend.compile(spec)  # `spec` and `weights` from the reservoir example above
print(compiled.dry_run_info())  # inspect offline, before submitting
handle = compiled.submit(weights, np.zeros(spec.n_inputs))
print(handle.status())
dist = handle.result()  # blocks until the cloud job completes

The access token is read from the token= argument, then from $QUANDELA_TOKEN or $PERCEVAL_TOKEN, and then from the token that Perceval itself would use (saved with pcvl.RemoteConfig). It is never written in the code. The processor connects once per compiled handle, and the connection is reused for later runs. run and result() wait at most poll_timeout_seconds (3600 by default). They raise an error if the job fails, or if the post-selection on the number of detected photons (min_detected_photons) discards every shot. See examples/quandela_cloud_hardware.py, and the live test marked qpu in tests/photonics/test_quandela_cloud.py, which you run with pytest -m qpu and a token set.

How it maps to the qubit core

The photonics package has the same layers as the core, for a different model of computation:

Qubit core	Photonics package
`core/ir.py`	`photonics/ir.py`
`core/provider.py` (the protocol, and the general parameter-shift rule used when a backend has nothing faster)	the `photonics` executor and `photonics/gradients.py`
`core/backends/*_backend.py`	the Perceval executor, which is a photonic backend
`torch/` (the PyTorch interface)	`PhotonicCircuit` (the photonic PyTorch interface)

The Perceval executor is a photonic backend, not a qubit Provider. The PyTorch adapter belongs to the photonics package, not to the qubit torch/ package. The Parameter tree, the separation between compiling and running, and the autograd.Function bridge are reused from the core without changes.

Status and limitations

Local simulation uses the strong linear optical simulation (SLOS) algorithm of Perceval. It is exact within the sizes it can handle. Its cost grows with the numbers of photons and modes, and the sizes it can handle are smaller than for a qubit statevector of size 2ⁿ.
Execution on Quandela Cloud, on a real QPU or a cloud simulator, goes through QuandelaCloudBackend. The results are sampled, so there are no gradients.
One scalar observable is evaluated at a time, for one sample. There is no batching yet.
Only observables that are diagonal in the Fock basis are supported.

Roadmap

These items are planned, in no fixed order:

A PhotonicMeasurement output type will allow several observables per evaluation and readouts that are not diagonal. It will also report the success probability of the post-selection with a sampled distribution.
run, expectation and PhotonicCircuit will accept a leading batch axis.
The Perceval circuit will be built once, with symbolic pcvl.P parameters that are bound at each call, in place of being rebuilt at every call.
Photonic reservoirs will get a quantum memory (an analogue of Carry), which needs Fock states to be injected between time steps.
QuandelaCloudBackend will move off RemoteProcessor and RemoteJob, which Perceval 1.3 keeps under perceval.runtime.legacy (they are still re-exported as pcvl.*).
MerLin, the PyTorch framework of Quandela for quantum machine learning, will be evaluated as an alternative to the parameter-shift bridge, with native autograd.