Hey! I am trying to sample from a graph’s covariance matrix from a GBS. However I keep getting that it is not a classical covariance matrix as output. Here is my code:-

A = nx.adjacency_matrix(gnp_random_graph(10,0.5)).todense()
_, s, _ = np.linalg.svd(A, full_matrices=True)
c = 1 / ( np.max(s) + 1e-8 )
Ab = c * np.block([
[A, np.zeros(A.shape)],
[np.zeros(A.shape), A]
])
X = np.block([
[np.zeros(A.shape), np.identity(A.shape[0])],
[np.identity(A.shape[0]), np.zeros(A.shape)]
])
I = np.identity(Ab.shape[0])
Q = np.linalg.inv(I - X @ Ab)
Q=Q.round(decimals=3)
cov=Q-I/2
print(cov)
thewalrus.samples.hafnian_sample_classical_state(cov,samples=1)

I have tried all sampling algorithms and none work, could you look into this?

Hey @MUSHKAN_SUREKA! Digging into the source code for thewalrus.samples.hafnian_sample_classical_state(), the function thewalrus.quantum.is_classical_cov() gets called and returns False from a call to thewalrus.quantum.is_valid_cov(). Looking into the various reasons why thewalrus.quantum.is_valid_cov() could return False, it appears that your matrix cov violates the the uncertainty relation!

I copy-pasted the source code from thewalrus.quantum.is_valid_cov() to make sure that this is indeed the problem:

I actually coded the covariance matrix wrt the Xanadu paper (https://arxiv.org/pdf/1905.12646.pdf) equation.2. I do not think there is an error in the code since I have used it for other purposes. Could you please look into it?

It looks like the convention in Maria’s paper is to subtract \mathbb{I}/2 from the covariance matrix Q. In thewalrus, this is baked-in! So, we just need to set cov = Q rather than cov = Q - 0.5 * I.

Covariance matrices need to be scaled by hbar. Typically, hbar = 1 is the convention. But, in our photonics libraries — mrmustard and thewalrus — the default value of hbar is set to 2 (usually hbar appears with a factor of 1/2, so hbar / 2 = 1 ). So you also need to scale cov: cov = hbar*Q. Having said that, if you’re like me and prefer working with hbar = 1, then just make sure you specify that whenever you call a thewalrus function.

Your code works now

import networkx as nx
import thewalrus
import numpy as np
hbar = 1
A = nx.adjacency_matrix(nx.gnp_random_graph(10, 0.5)).todense()
_, s, _ = np.linalg.svd(A, full_matrices=True)
c = 1 / (np.max(s) + 1e-8)
Ab = c * np.block([[A, np.zeros(A.shape)], [np.zeros(A.shape), A]])
X = np.block(
[
[np.zeros(A.shape), np.identity(A.shape[0])],
[np.identity(A.shape[0]), np.zeros(A.shape)],
]
)
I = np.identity(Ab.shape[0])
Q = np.linalg.inv(I - np.dot(X, Ab))
cov = hbar*Q
thewalrus.samples.hafnian_sample_classical_state(cov, samples=1, hbar=hbar)

I also removed Q.round(decimals=3) . Not sure why that’s needed .