Boson Sampling Applications with Strawberry Fields

Gaussian Boson Sampling (GBS) is more than just a demonstration of quantum computational advantage. The sf.apps module in Strawberry Fields provides tools that connect GBS to practical problems in chemistry and graph theory. In this tutorial, you will use GBS to simulate molecular vibronic spectra (computing Franck-Condon factors) and to measure graph similarity using graph kernels derived from photon counting statistics.

Prerequisites

Install Strawberry Fields along with The Walrus, Xanadu’s library for hafnian calculations:

pip install strawberryfields thewalrus numpy matplotlib

Why Gaussian Boson Sampling Is Classically Hard

In a GBS experiment, squeezed light is sent through a linear interferometer, and the output photon numbers are measured. The probability of observing a particular photon counting pattern is related to the hafnian of a submatrix of the interferometer’s transfer matrix. Computing hafnians is a #P-hard problem, meaning that no polynomial-time classical algorithm is known.

More precisely, for an $N$-mode GBS device with output pattern $\mathbf{n} = (n_1, n_2, \ldots, n_N)$, the probability is:

$P(\mathbf{n}) \propto \frac{|\text{Haf}(A_{\mathbf{n}})|^2}{\prod_i n_i!}$

where $A_{\mathbf{n}}$ is a submatrix of the adjacency matrix determined by the output pattern, and $\text{Haf}$ is the hafnian function. While computing this probability exactly is classically intractable, a photonic device naturally samples from this distribution.

Xanadu demonstrated this in 2022 with Borealis, a 216-mode photonic processor that sampled from a GBS distribution in seconds, a task estimated to require 9,000 years on the best classical supercomputers.

Application 1: Molecular Vibronic Spectra

Background

When a molecule absorbs a photon and undergoes an electronic transition, its vibrational state also changes. The resulting vibronic spectrum depends on the overlap between vibrational wavefunctions of the ground and excited electronic states. These overlaps are called Franck-Condon factors, and they determine the intensity of each spectral line.

Computing Franck-Condon factors is equivalent to computing the probability amplitudes of a specific bosonic transformation, which maps directly onto GBS. The Duschinsky rotation matrix (relating normal modes of the ground and excited states) becomes the interferometer unitary, and the displacement between equilibrium geometries becomes the squeezing and displacement parameters.

Setting Up a Molecular Simulation

The sf.apps.vibronic module provides functions to convert molecular parameters into a GBS circuit. You need three pieces of molecular data:

• The Duschinsky rotation matrix $U$ (relating normal modes between electronic states)
• The displacement vector $\delta$ (shift in equilibrium geometry)
• The frequency vectors $\omega$ and $\omega'$ (vibrational frequencies in ground and excited states)

For this tutorial, we will set up a simple two-mode model that captures the essential physics.

import numpy as np
from strawberryfields.apps import vibronic, sample

# Define molecular parameters for a simple 2-mode model
# Duschinsky rotation matrix (relates normal modes between states)
U_duschinsky = np.array([
    [0.9, 0.1],
    [-0.1, 0.9]
])
# Normalize to make it unitary (approximate for illustration)
from numpy.linalg import svd
u, s, vh = svd(U_duschinsky)
U_duschinsky = u @ vh

# Displacement vector (shift in equilibrium geometry)
delta = np.array([0.5, 0.3])

# Vibrational frequencies (in wavenumbers, cm^-1)
omega_ground = np.array([1000.0, 2000.0])   # Ground state frequencies
omega_excited = np.array([1100.0, 1900.0])   # Excited state frequencies

# Temperature (Kelvin) affects the initial thermal population
T = 300.0  # Room temperature

# Convert molecular parameters to GBS parameters
gbs_params = vibronic.gbs_params(
    w=omega_ground,
    wp=omega_excited,
    Ud=U_duschinsky,
    delta=delta,
    T=T
)

# gbs_params contains: t (two-mode squeezing), U1, r (local squeezing),
# U2 (interferometer), alpha (displacement)
t, U1, r, U2, alpha = gbs_params

print("Two-mode squeezing parameters:", t)
print("Local squeezing parameters:", r)
print("Displacement parameters:", alpha)
print("Interferometer U1 shape:", U1.shape)
print("Interferometer U2 shape:", U2.shape)

Sampling Franck-Condon Factors

Now we generate samples from the GBS circuit. Each sample corresponds to a particular vibrational transition, and the frequency of each sample approximates the corresponding Franck-Condon factor.

import strawberryfields as sf
from strawberryfields.ops import (
    S2gate, Interferometer, Sgate, Dgate, MeasureFock
)

num_modes = len(omega_ground)

# Build the GBS circuit from molecular parameters
prog = sf.Program(num_modes * 2)  # Need 2N modes for the vibronic circuit

with prog.context as q:
    # Two-mode squeezing between mode pairs
    for i in range(num_modes):
        S2gate(t[i]) | (q[i], q[i + num_modes])

    # First interferometer on the first N modes
    Interferometer(U1) | q[:num_modes]

    # Local squeezing on first N modes
    for i in range(num_modes):
        Sgate(r[i]) | q[i]

    # Second interferometer on the first N modes
    Interferometer(U2) | q[:num_modes]

    # Displacement on first N modes
    for i in range(num_modes):
        Dgate(np.abs(alpha[i]), np.angle(alpha[i])) | q[i]

    # Measure photon numbers in all modes
    MeasureFock() | q

eng = sf.Engine("fock", backend_options={"cutoff_dim": 6})
result = eng.run(prog, shots=10000)

samples = result.samples
print(f"Generated {len(samples)} samples")
print(f"Sample shape: {samples.shape}")
print(f"First 5 samples:\n{samples[:5]}")

Computing the Vibronic Spectrum

Each sample represents a vibrational transition. The energy of the transition is determined by the photon numbers in each mode, and the frequency of occurrence approximates the Franck-Condon factor for that transition.

import matplotlib.pyplot as plt

# Compute transition energies from samples
# Energy = sum of (n_i * omega_i) for excited state modes
# minus sum of (m_i * omega_i) for ground state modes (thermal)

# For simplicity, compute the excited state vibrational energy
energies = []
for s in samples:
    # First N modes correspond to the excited state
    excited_quanta = s[:num_modes]
    # Last N modes correspond to the ground state
    ground_quanta = s[num_modes:]

    # Transition energy (in cm^-1)
    e_excited = np.sum(excited_quanta * omega_excited)
    e_ground = np.sum(ground_quanta * omega_ground)
    energies.append(e_excited - e_ground)

energies = np.array(energies)

# Plot the vibronic spectrum
plt.figure(figsize=(10, 5))
plt.hist(energies, bins=50, density=True, alpha=0.7, color="steelblue",
         edgecolor="black")
plt.xlabel("Transition Energy (cm$^{-1}$)")
plt.ylabel("Intensity (normalized)")
plt.title("Simulated Vibronic Spectrum via GBS")
plt.grid(True, alpha=0.3)
plt.savefig("vibronic_spectrum.png", dpi=150)
plt.show()

# Print the most common transitions
from collections import Counter
transition_counts = Counter(map(tuple, samples))
most_common = transition_counts.most_common(10)
print("\nMost common transitions (photon pattern: count):")
for pattern, count in most_common:
    print(f"  {pattern}: {count} ({100*count/len(samples):.1f}%)")

The resulting histogram approximates the vibronic absorption spectrum of the molecule. Peaks in the histogram correspond to strong Franck-Condon factors, meaning those vibrational transitions have high probability.

Application 2: Graph Similarity via GBS

Background

Graphs appear everywhere in science: molecular structures, social networks, communication topologies. Measuring how similar two graphs are is a fundamental problem. Classical graph kernels exist but often struggle with large graphs. GBS provides a novel approach to graph similarity.

The connection works as follows. Given a graph with adjacency matrix $A$, you can program a GBS device so that the sampling probabilities encode properties of the graph. Specifically, the probability of detecting photons in a subset of modes relates to the hafnian of the corresponding submatrix of $A$, which counts the number of perfect matchings in that subgraph. By comparing the photon counting statistics of two GBS devices programmed with different graphs, you obtain a similarity measure.

Encoding Graphs for GBS

The sf.apps.graph module (in older versions) and sf.apps.similarity module provide tools for this. Here we will build the graph encoding manually to illustrate the underlying mechanism.

import numpy as np
import strawberryfields as sf
from strawberryfields.ops import Sgate, Interferometer, MeasureFock

def adjacency_to_gbs(A, mean_photon=2.0):
    """
    Convert a graph adjacency matrix to GBS circuit parameters.

    The GBS device is programmed so that the squeezing parameters
    and interferometer encode the graph structure.
    """
    from numpy.linalg import svd, eig

    # Eigendecomposition of the adjacency matrix
    eigenvalues, eigenvectors = eig(A)
    eigenvalues = np.real(eigenvalues)
    eigenvectors = np.real(eigenvectors)

    # Set squeezing parameters proportional to eigenvalues
    # Scale so that mean photon number matches target
    max_eig = np.max(np.abs(eigenvalues))
    if max_eig > 0:
        scale = np.tanh(mean_photon / len(A)) / max_eig
        squeezing = np.arctanh(np.clip(np.abs(eigenvalues) * scale, 0, 0.99))
    else:
        squeezing = np.zeros(len(A))

    # Signs of eigenvalues determine squeezing phases
    phases = np.where(eigenvalues >= 0, 0.0, np.pi)

    return eigenvectors, squeezing, phases


def sample_gbs_from_graph(A, num_samples=5000, mean_photon=2.0, cutoff=6):
    """
    Run GBS sampling for a graph with adjacency matrix A.
    """
    num_modes = len(A)
    U, sq_r, sq_phi = adjacency_to_gbs(A, mean_photon)

    prog = sf.Program(num_modes)
    with prog.context as q:
        # Apply squeezing encoding the graph eigenvalues
        for i in range(num_modes):
            Sgate(sq_r[i], sq_phi[i]) | q[i]

        # Apply interferometer encoding the graph eigenvectors
        Interferometer(U) | q

        # Measure photon numbers
        MeasureFock() | q

    eng = sf.Engine("fock", backend_options={"cutoff_dim": cutoff})
    result = eng.run(prog, shots=num_samples)
    return result.samples

Comparing Two Graphs

Let us define two graphs and compare them using their GBS sampling distributions.

import matplotlib.pyplot as plt
from collections import Counter

# Graph 1: A cycle graph on 4 nodes (square)
A1 = np.array([
    [0, 1, 0, 1],
    [1, 0, 1, 0],
    [0, 1, 0, 1],
    [1, 0, 1, 0]
], dtype=float)

# Graph 2: A path graph on 4 nodes (line)
A2 = np.array([
    [0, 1, 0, 0],
    [1, 0, 1, 0],
    [0, 1, 0, 1],
    [0, 0, 1, 0]
], dtype=float)

# Graph 3: A complete graph on 4 nodes
A3 = np.array([
    [0, 1, 1, 1],
    [1, 0, 1, 1],
    [1, 1, 0, 1],
    [1, 1, 1, 0]
], dtype=float)

# Sample from each graph
print("Sampling from Graph 1 (cycle)...")
samples1 = sample_gbs_from_graph(A1, num_samples=5000)
print("Sampling from Graph 2 (path)...")
samples2 = sample_gbs_from_graph(A2, num_samples=5000)
print("Sampling from Graph 3 (complete)...")
samples3 = sample_gbs_from_graph(A3, num_samples=5000)

Computing the Graph Kernel

The graph kernel compares the photon counting distributions. We compute a feature vector from the GBS samples, where each element corresponds to the probability of a particular photon number orbit (the total photon count, ignoring which modes detected them).

def photon_number_distribution(samples, max_photons=8):
    """
    Compute the distribution over total photon numbers from GBS samples.
    This serves as a feature vector for the graph.
    """
    total_photons = np.sum(samples, axis=1)
    distribution = np.zeros(max_photons + 1)
    for n in range(max_photons + 1):
        distribution[n] = np.mean(total_photons == n)
    return distribution


def graph_kernel(dist1, dist2):
    """
    Compute the similarity between two graphs using their
    photon number distributions. Uses the fidelity (Bhattacharyya
    coefficient) as the kernel.
    """
    return np.sum(np.sqrt(dist1 * dist2))


# Compute feature vectors
dist1 = photon_number_distribution(samples1)
dist2 = photon_number_distribution(samples2)
dist3 = photon_number_distribution(samples3)

# Compute pairwise similarities
k12 = graph_kernel(dist1, dist2)
k13 = graph_kernel(dist1, dist3)
k23 = graph_kernel(dist2, dist3)

print(f"Kernel(cycle, path):     {k12:.4f}")
print(f"Kernel(cycle, complete): {k13:.4f}")
print(f"Kernel(path, complete):  {k23:.4f}")

# Visualize the photon number distributions
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
names = ["Cycle (C4)", "Path (P4)", "Complete (K4)"]
dists = [dist1, dist2, dist3]

for ax, name, dist in zip(axes, names, dists):
    ax.bar(range(len(dist)), dist, color="steelblue", edgecolor="black")
    ax.set_xlabel("Total Photon Number")
    ax.set_ylabel("Probability")
    ax.set_title(name)
    ax.set_ylim(0, max(d.max() for d in dists) * 1.2)

plt.tight_layout()
plt.savefig("graph_gbs_distributions.png", dpi=150)
plt.show()

The kernel values reveal structural relationships between the graphs. The cycle and path graphs (both sparse, both regular or near-regular) tend to have more similar distributions than either has with the complete graph. This reflects the fact that the GBS sampling distribution encodes structural properties of the graph, such as the number and size of perfect matchings.

Using sf.apps for Graph Problems

Strawberry Fields provides higher-level functions in the sf.apps module that automate much of this workflow. The sf.apps.sample module generates GBS samples from an adjacency matrix, and sf.apps.similarity computes graph kernels directly.

from strawberryfields.apps import sample as gbs_sample

# Generate GBS samples directly from an adjacency matrix
# This handles all the squeezing and interferometer setup internally
samples_auto = gbs_sample.sample(A1, n_mean=2.0, n_samples=5000)
print(f"Auto-generated {len(samples_auto)} samples from graph A1")
print(f"Sample: {samples_auto[0]}")

The sf.apps.similarity module provides ready-made kernel functions:

from strawberryfields.apps import similarity

# Compute feature vectors based on different orbit types
# An "orbit" groups photon patterns by their total photon count
feature1 = similarity.feature_vector_sampling(samples_auto, event_photon_numbers=[1, 2, 3, 4])
print(f"Feature vector for graph 1: {feature1}")

Practical Considerations

Sampling noise. Because GBS is a sampling-based method, all computed quantities carry statistical uncertainty. More samples yield better estimates, but the number of samples needed grows with the precision required. For the graph kernel application, a few thousand samples typically suffice for small graphs, but larger graphs may require significantly more.

Classical simulation limits. Classical simulation of GBS via hafnian computation is feasible for up to about 40 to 50 modes. Beyond that, only approximate methods or actual photonic hardware can produce samples efficiently. The thewalrus library provides optimized hafnian computations for classical validation of small instances.

Fock space cutoff. When using the fock backend, the cutoff_dim parameter truncates the infinite-dimensional Fock space. If too small, the simulation will be inaccurate because high-photon-number states are cut off. As a rule of thumb, set the cutoff to at least twice the expected mean photon number per mode.

Hardware access. For problems requiring more than a few modes, Xanadu’s cloud-accessible Borealis processor offers 216 squeezed modes. Access is available through the Xanadu Cloud platform, and Strawberry Fields provides the sf.RemoteEngine interface for submitting jobs.

Summary

Gaussian Boson Sampling connects the abstract problem of sampling from photonic circuits to concrete applications in chemistry and graph theory. The vibronic spectra application exploits the mathematical equivalence between Franck-Condon factors and GBS probabilities, allowing photonic processors to simulate molecular spectroscopy. The graph similarity application encodes graph structure into the GBS device and extracts structural features from the photon counting statistics. Both applications are available through the sf.apps module, which provides high-level interfaces that abstract away the circuit construction details. As photonic hardware scales beyond what classical computers can simulate, these applications represent some of the most promising paths toward practical quantum advantage.