Feature request for ThresholdBosonSampling

[GongSiqiu]

Your library is excellent and the sampling speed is much faster than thewalrus! It would be better if it had the following features.

1. Set the upbound for sampling. It is known that it will take a lot of time if the simulation involves too many photons. Sometimes, it is only needed to generate samples whose photon-click-number are less than the upbound. It will save much time if it can skip the calculation for photons more than the upbound.
2. Use a new algorithm. There is a better algorithm for generating GBS samples with threshold detectors. It will have an advantage if there is a high probability of collisions between photons. In addition, the algorithm significantly reduces the computation complexity to calculate Hafnian with several photons in one port.

[GregoryMorse]

Hello, as for the new repeated Eigenvalue trace algorithm, we already have discussed it and have plans for it. In fact, the code needs some generalization and is a big displaced so its for GBS but not for raw (loop) Hafnian calculations.

So we need to integrate this code: https://github.com/Budapest-Quantum-Computing-Group/piquassoboost/blob/main/piquassoboost/sampling/simulation_strategies/source/RepeatedColumnPairs.cpp

Into: https://github.com/Budapest-Quantum-Computing-Group/piquassoboost/blob/main/piquassoboost/sampling/source/PowerTraceHafnianRecursive.cpp and https://github.com/Budapest-Quantum-Computing-Group/piquassoboost/blob/main/piquassoboost/sampling/source/PowerTraceLoopHafnianRecursive.cpp

While generalizing the occupancy to a per row/column not per pair basis.

But fundamentally we still use binomial coefficients and the Eigenvalue trace method AFAICT.

I will be revisiting this particular issue and slightly improving the accuracy on all Hafnian calculators soon. Pinging @rakytap on this.

[GongSiqiu]

I add a parameter "upbond" for ThresholdBosonSampling::simulate. It does save me a lot of time. However, I can just hard code it and can't set it in python. What should I do to manage it?

[rakytap]

Dear GongSiqiu, Please push up a new branch containing your modification and a python script showing us the way you are making the calls, we will expand the python interface for your need.

[GongSiqiu]

I don't know how to push up a new branch to your repertory. These are the modified codes and expected python call upper-bounds.zip

[rakytap]

Dear GongSiqiu,

You can access the modiofied code at branch rp_upbound. Since we dont support in piquasso the overall upper bound on the photon count, you need to access this functionality through low-level API:

from piquassoboost.sampling.simulation_strategies import ThresholdBosonSampling

covariance_mtx = np.real( ... some covariance matrix ... ) # needs to be real

th = ThresholdBosonSampling.ThresholdBosonSampling( interferometer_param )

th.set_photon_upper_bound(30)

result = th.simulate(shots)

Piquaaso also offers the ability to chose between hafnian based or torontonian based simulation of threshold BS. In the hafnian variant we still need to finish some small development work to account for photon multiplicities, but it is competitive with the torontonian based one, especially at higher photon counts.

Here is a simple example to chose between the hafnian/torontonian based variant (see line )

import numpy as np import piquasso as pq import piquassoboost as pqb from piquassoboost.config import BoostConfig import thewalrus cov = np.load("./cov.npy")

d = 50 shots = 10000

config = BoostConfig(use_torontonian=True)

with pq.Program() as gaussian_boson_sampling: pq.Q(all) | pq.Covariance(thewalrus.symplectic.xxpp_to_xpxp(cov)) pq.Q(all) | pq.ThresholdMeasurement()

simulator = pqb.BoostedGaussianSimulator(d=d, config=config)

result = simulator.execute(gaussian_boson_sampling, shots=shots) samples = np.array(result.samples)

Dont forget to recompile the code after checking out the rp_upbound branch

[rakytap]

Also, today we updated the compatibility with the newer version of piquasso, so probably you will need to update also the piquasso submodule by 'git submodule update'

[GongSiqiu]

Do I need to modify the code at branch rp_upbound? It seems that this branch is not changed recently. I manage to compile it and try your python code. However, it tells me that 'ThresholdBosonSampling' object has no attribute 'set_photon_upper_bound'

[rakytap]

Yeah, I forgot to push up the git commit. Please update your local repository.

[GongSiqiu]

The sampling result seems wrong. Did I make some mistakes?

from piquassoboost.sampling.simulation_strategies import ThresholdBosonSampling

covariance_mtx = np.real(thewalrus.symplectic.xxpp_to_xpxp(cov)) # needs to be real

th = ThresholdBosonSampling.ThresholdBosonSampling(covariance_mtx)
shots = 1000

th.set_photon_upper_bound(10)

samples = np.array(th.simulate(shots))

[GongSiqiu]

The click-number distribution looks like this

[rakytap]

Would not

samples = np.array(th.simulate(shots)) plt.plot( samples.sum( axis = 1) )

make things right?

[GongSiqiu]

I noticed that mistake so I deleted that reply :) That is not the major problem. The TVD between samples and theory is too big.

[rakytap]

Well, I think we would need to know more details about your work in order to help you. (You know, I only added a python interface call to set the upbound, while taking your modifications of the sampling code with no modifications.)

[GongSiqiu]

I just randomly generate a 20*20 covariance matrix which is the result of feeding single mode squeezed vaccum states to a unitary interferometer and do the sampling.

Is there a default upbound? If I generate samples in the original way, the result is correct while it would be wrong if I use the low-level API.

[rakytap]

1) ok, we need to define what does the "original way" mean. We support two variants of threshold BS: 1) hafnian based, 2) torontonian based. It is possible to switch between the two engines at high level API by overriding the default settings via config = BoostConfig(use_torontonian=True/False) (see my comment above)

So the question is which one did you use in your original experiments?

2) Did you verify your hardcoded modification in the source code?

3) a difference between the results generated by low and high level APIs might be coming from the input matrix used to initialize the C++ ThresholdBosonSampling class. This should be checked. If these input matrices are identical, then the low/high API calls should give identical results.

[rakytap]

Or put it another way: if you provide us with two example scripts that should give the same results then it is possible to start debugging.

[GongSiqiu]

This is the "original way"

import piquasso as pq
import piquassoboost as pqb
from piquassoboost.config import BoostConfig
import thewalrus
cov = G.cov

d = 10
shots = 100000

# config = BoostConfig(use_torontonian=True)
config = BoostConfig(use_torontonian=True)

with pq.Program() as gaussian_boson_sampling:
    pq.Q(all) | pq.Covariance(thewalrus.symplectic.xxpp_to_xpxp(cov))
    pq.Q(all) | pq.ThresholdMeasurement()

simulator = pqb.BoostedGaussianSimulator(d=d, config=config)

result = simulator.execute(gaussian_boson_sampling, shots=shots)
samplesOri = np.array(result.samples)

In this way, the hardcoded upbound do take effect.

These input matrices are identical

[GongSiqiu]

debug.zip This an example

[rakytap]

At first glance I see two differences: 1) In the low level API the xxpp variant of the covariance matrix should be used (see https://github.com/Budapest-Quantum-Computing-Group/piquassoboost/blob/main/piquassoboost/gaussian/calculations.py#L66)

2) the covariance matrix should be scaled by 2*hbar (see https://github.com/Budapest-Quantum-Computing-Group/piquassoboost/blob/main/piquassoboost/gaussian/calculations.py#L67)

By default hbar=1, I think.

[GongSiqiu]

It works! Thank you very much!! I have tried other hbar and xp order combinations but lost this pair.

[rakytap]

You are welcome!