Fock backend of localsimulator and thewalrus decomposition not giving the desirable result

[minseok1999]

Before posting a bug report

• I have searched exisisting GitHub issues to make sure the issue does not already exist.

Expected behavior

I wanted to follow the full steps of encoding an adjacency matrix of a graph into an interferometric setup and single mode squeezed states for Gaussian boson sampling experiment. To do this, I first used thewalrus library to decompose the adjacency matrix using Takagi decomposition. After I have decomposed the matrix, I have then used rectangular_phase_end from the strawberry fields library to convert my unitary gate into a set of beam splitters. Finally, I have used the Fock backend of the local simulator engine and fed in the squeezing parameters and Beam splitter setup I have generated to get the photo count result corresponding to my specific adjacency matrix.

I expected the photocount pattern of [111100] to appear for the largest number of times since that pattern corresponds to the 4 dimensional subgraph choice of maximum Hafnium value to my 6 dimensional adjacency matrix. (Hafnium value of 3, Hafnium squared value of 9 whereas other choices such as [110011] only has single perfect matching) I have tested my graph using hafnium_sample_graph from thewalrus library as well, and I got 231 samples of [111100] among 10000 samples, which is quite desirable.

Actual behavior

What I got after postselecting only those with ones and zeros (subspace used for the estimation of hafnian) was something like this,
[1 1 1 0 0 1]
[1 0 0 1 1 1]
[1 0 1 0 1 1]
[0 1 1 1 0 1]
[1 1 0 1 1 0]
[1 1 1 0 1 0]
[0 0 1 1 1 1]
[1 1 1 0 0 1]

only 8 cases out of 1000 samples and I could not see the answer case,i.e, [111100] ,to the maximum hafnian problem at all.

What's more, mean photon number at the output was around 0.97 whereas I aimed it at 4 to best see the 4 dimensional subgraph among my 6 dimensional adjacency matrix.

Reproduces how often

Reproduces nearly every time.

System information

Strawberry Fields: a Python library for continuous-variable quantum circuits.
Copyright 2018-2020 Xanadu Quantum Technologies Inc.

Python version:            3.10.9
Platform info:             macOS-13.4.1-arm64-arm-64bit
Installation path:         /Users/ryuminseok/anaconda3/lib/python3.10/site-packages/strawberryfields
Strawberry Fields version: 0.23.0
Numpy version:             1.23.5
Scipy version:             1.10.0
SymPy version:             1.12
NetworkX version:          3.1
The Walrus version:        0.20.0
Blackbird version:         0.5.0
XCC version:               0.3.0
TensorFlow version:        None

Source code

import numpy as np

matrix = np.array(
[[1,1,1,1,0,0],
[1,1,1,1,0,0],
[1,1,1,1,0,1],
[1,1,1,1,0,0],
[0,0,0,0,1,1],
[0,0,1,0,1,1]]
)-np.eye(6)
print(matrix)

import matplotlib.pyplot as plt
import networkx as nx
# Create the entire graph
gr = nx.Graph(matrix)

default_edge_color = 'gray'

edge_colors = {edge: default_edge_color for edge in gr.edges()}

node_color = 'orange'
nx.draw(gr, node_size=500,node_color=node_color, edge_color=[edge_colors[edge] for edge in gr.edges()],with_labels=True)

# Display the graph
plt.show()

import strawberryfields as sf

dim=np.shape(matrix)
print(dim[0])
"write the mean photon number you want to give at the output"
meanphoton=4
meanphotonpermode=meanphoton/dim[0]
print(meanphotonpermode)
result=sf.decompositions.graph_embed(matrix,meanphotonpermode)
squeeze=result[0]
uni=result[1]

setup=sf.decompositions.rectangular_phase_end(uni)
newt=setup[0]
"Print out the beam splitters"
print(np.round(newt,2))
ss=np.shape(newt)[0]

from strawberryfields import ops

nr_modes=dim[0]
prog = sf.Program(nr_modes)

with prog.context as q:
    for i in range(dim[0]) :
        ops.Sgate(squeeze[i]) |(q[i])
    for j in range(ss) :
        ops.BSgate(np.round(newt[j][2],2),np.round(newt[j][3],2)) |(q[newt[j][0]],q[newt[j][1]])

    ops.MeasureFock() | q


from strawberryfields import LocalEngine

eng = sf.Engine("fock", backend_options={"cutoff_dim": 4})
"Check if circuit is working as desired"
eng.run(prog)
eng.print_applied()

from pytictoc import TicToc

t=TicToc()
t.tic()

sample_array = []  # Array to store the samples
"Type in How many samples you want to put"
num=1000
for _ in range(num):
    results = eng.run(prog)
    sample=results.samples
    sample_array.append(sample)  # Add the results to the array

t.toc()

sample_array = np.array(sample_array)

peeled_array = np.array(sample_array)[:, 0]

mean_count = np.mean(peeled_array.sum(axis=1))
"What is the mean photon number of the entire output?"
print(mean_count)

concatenated_array = np.concatenate(sample_array, axis=0)
"Postselect only those with 0 or 1"
filtered_array = [arr for arr in concatenated_array if all(elem in {0, 1} for elem in arr)]
"Postselect 4 dim subgraph"
postselected_array = [arr for arr in filtered_array if np.sum(arr) == 4]

for arr in postselected_array:
    print(arr)

"Count the number of answer samples to MaxHaf problem"
target_sequence = np.array([1, 1, 1, 1, 0, 0])
count = 0

for arr in postselected_array:
    if np.array_equal(arr, target_sequence):
        count += 1

print(count) 

Tracebacks

No response

Additional information

I add information of Hafnian value of every 4 dimensional subgraph to my specific graph
 
Submatrix:

i)Perfect matching

ii) (prob)∝|HafA_s |^2


0123  (a case)

i) 3

ii) 9

xxx5  (10 cases)

i) 1

ii) 1

xxx4  (4 cases)

i) 0

ii) 0

[isaacdevlugt]

Hey @minseok1999! Apologies for the delayed response. Can you distill your code down to something minimal that still reproduces the behaviour? It would help me greatly 😄

[minseok1999]

Hey @minseok1999! Apologies for the delayed response. Can you distill your code down to something minimal that still reproduces the behaviour? It would help me greatly 😄

Thank you for considering my issue again! I have edited with the distilled source code as you requested.

[isaacdevlugt]

@minseok1999 can you surround your code with three backticks? It will render it like this, which github will let me copy-paste.

# Code here

You can also get rid of the notebook block stuff (e.g., In[1]).

[minseok1999]

@minseok1999 can you surround your code with three backticks? It will render it like this, which github will let me copy-paste.

# Code here

You can also get rid of the notebook block stuff (e.g., In[1]).

I have removed notebook blocks and wrapped my code using backticks as you requested!

[isaacdevlugt]

Awesome thanks!

It might be something to do with eng = sf.Engine("fock", backend_options={"cutoff_dim": 4}) — if the cutoff dimension is 4, then each Fock representation of the components has a cutoff of 4 (i.e., it supports only up to 3 photons). If throughout the process a mode needs to populate more than 3 photons, it can’t. Let us know if that helps!

[minseok1999]

Awesome thanks!

It might be something to do with eng = sf.Engine("fock", backend_options={"cutoff_dim": 4}) — if the cutoff dimension is 4, then each Fock representation of the components has a cutoff of 4 (i.e., it supports only up to 3 photons). If throughout the process a mode needs to populate more than 3 photons, it can’t. Let us know if that helps!

I have raised the cutoff dimension up to 5 and 6 and it clearly gave higher photon 'mean_count' (average photon number at the output port out of 1000 samples) as 1.46 and 1.54. This is better but clearly not good enough cause I have aimed my mean photon number at 4 when I decomposed my graph into experimental setup of beam splitters and single mode squeezed states. But the problem is, I cannot raise this cutoff dimension above 7 because of memory deficiency of my MacBook Air. Is there a way to reliably simulate 6 mode Gaussian Boson Sampling experiment without this memory issue other than the local Fock basis simulator? I know hafnian_sample_graph gives gaussian boson sampling pattern but what I want is to simulate this from the experimental setup of beam splitters so that I can simulate in the future how this sampling pattern is affected upon attaching some other components.

[isaacdevlugt]

Is there a way to reliably simulate 6 mode Gaussian Boson Sampling experiment ...

These demos are worth checking out:

• https://strawberryfields.ai/photonics/apps/run_tutorial_sample.html
• https://strawberryfields.ai/photonics/apps/run_tutorial_training.html

Let us know if those help you get anywhere!

[minseok1999]

Is there a way to reliably simulate 6 mode Gaussian Boson Sampling experiment ...

These demos are worth checking out:

• https://strawberryfields.ai/photonics/apps/run_tutorial_sample.html
• https://strawberryfields.ai/photonics/apps/run_tutorial_training.html

Let us know if those help you get anywhere!

The first link you provided helped me with postselection stage after getting all the samples(tidier code than the one I wrote myself) but it didn't have much to do with full path optical simulation of Gaussian boson sampling(simulation done by writing out components manually and can be varied by attaching components myself) I am intending. And for the second one, I didn't scrutinize all the details yet but it looked like it was about some more sophisticated Gaussian boson sampling scheme of training the machine real time, which I don't need much at this stage.

I appreciate your advice though! If there is no such reliable simulator provided that is fit for me then I think I should build one myself or consider buying some gpu for more local kernel memory🙂

[isaacdevlugt]

It would definitely be worth it to see if raising the cutoff higher and higher helps you achieve the desired result. We can leave this issue open until the solution is determined 😄