Instability in quantum.probabilities about thewalrus HOT 9 CLOSED

Thanks @Marsll! We will take a look into this

from thewalrus.

Hi @Marsll --- Thanks for finding this. It is indeed a very mysterious bug. If you scan the phase of your coherent displacement and plot the normalization you get this

angles = np.arange(0,np.pi,0.01)
norm = [np.sum(probabilities(np.sqrt(2) * np.array([(alpha_abs * np.exp(1j * theta)).real, (alpha_abs * np.exp(1j * theta)).imag]) , cov, cutoff, hbar=1.0)) for theta in angles]
plt.plot(angles, norm)
plt.xlabel("displacement phase")
plt.ylabel("normalization")

On the other hand if the amplitude of your displacement is reduced to be |\alpha| = \sqrt{40} then things more or less work

I am tempted to say that it could be a numerical instability as you suggest, but then somehow it seems to orderly when you plot as a function of the angle.

One possible way to get around this issue is to use long double instead of double in the low level C++ functions that actually do the calculation.

from thewalrus.

Hi Nicolás, thanks so much for looking into this. It looks mysterious indeed.
Do I have any control over these C++ low level functions that you mention directly via the TheWalrus library, at the moment?

from thewalrus.

Hi Marcel --- You do. The chain of functions calls goes as follows

1. quantum.probabilities
2. quantum.state_vector
3. hafnian_batched (in _hermite_multidimensional.py)
4. hermite_multidimensional (in _hermite_multidimensional.py)
5. renorm_hermite_multidimensional

This last function lives in libwalrus.pyx where we Cython to be able to call renorm_hermite_multidimensional_cpp from hermite_multidimensional.hpp.

One worrisome aspect of your calculation is that ultimately all of this depends on a recursion relation starting from the vacuum amplitude of the Gaussian state. For the Gaussian state you are considering this quantity is (-1.1483424815398375e-15+3.522084812154078e-15j), i.e., dangerously close to machine epsilon.

from thewalrus.

Another thing I tried to test was if there was a communication problem between quantum.probabilities and the C++ code. I calculated <n|D(\alpha) S(r) |0> by first finding the Fock matrix representation of D(\alpha) S(r) which you can do as follows (using the same variables of your example code):

from thewalrus.quantum import fock_tensor
T = fock_tensor(S, np.array([alpha]), cutoff=100)
I can now plot this versus the result of probabilities as follows
plt.plot(np.abs(T[:,0])**2,"*")
plt.plot(probs_displaced_squeezed)

and I get

They agree perfectly and are both wrong (at least by a scale).

from thewalrus.

Finally, although the normalization is definitely wrong the probabilities are correct up to a scale. If you renormalize them
renorm_p = probs_displaced_squeezed / np.sum(probs_displaced_squeezed)
and used them to calculate the mean photon number and the variance
mean = (renorm_p) @ np.arange(cutoff)
var = var = (renorm_p) @ (np.arange(cutoff))**2 - mean**2
to get
print(mean,var)
50.271540042961355 36.29493009537737
you can easily verify that these are correct by running
photon_number_mean(mu, cov, 0, hbar=1)
and
photon_number_covar(mu, cov, 0,0, hbar=1)

which makes this even more mysterious....

from thewalrus.

Hi @nquesada any new updates/insights on this issue?

from thewalrus.

Yup, the problem is related to this line:

For @Marsll pref = (-1.1483424815398375e-15+3.522084812154078e-15j) while np.real_if_close(pref) = (-1.1483424815398375e-15+0j) . If the np.real_if_close is removed then the normalization comes out to be 0.9999999946493981 ~ 1 .

from thewalrus.

Solved in #215

from thewalrus.