Hi @CatalinaAlbornoz ,
I am more than happy to provide the full code of the _decompose _method, where I added a small code snippet! To make sure you can excatly reproduce what I encountered I will try to give you a step by step instructions what I did so that you can understand it.
- Start by running the following code:
import strawberryfields as sf
from strawberryfields.ops import *
import numpy as np
M = 4
# Create a random complex matrix
random_matrix = np.random.rand(M, M) + 1j * np.random.rand(M, M)
# Perform QR decomposition to get a unitary matrix
q, r = np.linalg.qr(random_matrix)
# Normalize to ensure unitary property (q should already be unitary, but we'll double-check)
unitary_matrix = q / np.linalg.norm(q, axis=0)
boson_sampling_triangular = sf.Program(4)
with boson_sampling_triangular.context as q:
# prepare the input fock states
Fock(1) | q[0]
Fock(1) | q[1]
Vac | q[2]
Fock(1) | q[3]
Interferometer(unitary_matrix, mesh = 'triangular', tol = 1e-10) | q
boson_sampling_triangular.compile(compiler="fock")
eng = sf.Engine(backend="fock", backend_options={"cutoff_dim": 5})
results = eng.run(boson_sampling_triangular)
# extract the joint Fock probabilities
probs_triangular = results.state.all_fock_probs()
#now try with mesh = 'rectangular'
boson_sampling_rectangular = sf.Program(4)
with boson_sampling_rectangular.context as q:
# prepare the input fock states
Fock(1) | q[0]
Fock(1) | q[1]
Vac | q[2]
Fock(1) | q[3]
Interferometer(unitary_matrix, mesh = 'rectangular', tol = 1e-10) | q
boson_sampling_rectangular.compile(compiler="fock")
eng = sf.Engine(backend="fock", backend_options={"cutoff_dim": 5})
results = eng.run(boson_sampling_rectangular)
# extract the joint Fock probabilities
probs_rect = results.state.all_fock_probs()
print(f"Probability to measure |1 1 0 1> for triangular mesh: {probs_triangular[1,1,0,1]}")
print(f"Probability to measure |1 1 0 1> for triangular mesh: {probs_rect[1,1,0,1]}")
print(f"{np.allclose(probs_triangular, probs_rect)}")
Here the results will differ between the rectangular and the triangular decomposition. But this should not be, as we still want to apply the same unitary matrix.
Now my fix for this problem:
- Change the
_decompose method to the following:
def _decompose(self, reg, **kwargs):
# pylint: disable=too-many-branches
mesh = kwargs.get("mesh", self.mesh)
tol = kwargs.get("tol", self.tol)
drop_identity = kwargs.get("drop_identity", self.drop_identity)
cmds = []
if mesh == "rectangular_compact":
phases = dec.rectangular_compact(self.p[0], rtol=tol, atol=tol)
cmds = _rectangular_compact_cmds(reg, phases)
elif mesh == "triangular_compact":
phases = dec.triangular_compact(self.p[0], rtol=tol, atol=tol)
cmds = _triangular_compact_cmds(reg, phases)
elif mesh == "sun_compact":
parameters, global_phase = dec.sun_compact(self.p[0], rtol=tol, atol=tol)
cmds = _sun_compact_cmds(reg, parameters, global_phase)
elif mesh == "triangular":
decomp_fn = getattr(dec, mesh)
BS1, R, BS2 = decomp_fn(self.p[0], tol=tol)
for n, expphi in enumerate(R):
# local phase shifts
q = np.log(expphi).imag if np.abs(expphi - 1) >= _decomposition_tol else 0
if not (drop_identity and q == 0):
cmds.append(Command(Rgate(np.mod(q, 2 * np.pi)), reg[n]))
for n, m, theta, phi, _ in BS1:
theta = theta if np.abs(theta) >= _decomposition_tol else 0
phi = phi if np.abs(phi) >= _decomposition_tol else 0
if not (drop_identity and theta == 0):
cmds.append(Command(BSgate(-theta, 0), (reg[n], reg[m])))
# Clements style beamsplitters
if not (drop_identity and phi == 0):
cmds.append(Command(Rgate(-phi), reg[n]))
elif not self.identity or not drop_identity:
decomp_fn = getattr(dec, mesh)
BS1, R, BS2 = decomp_fn(self.p[0], tol=tol)
for n, m, theta, phi, _ in BS1:
theta = theta if np.abs(theta) >= _decomposition_tol else 0
phi = phi if np.abs(phi) >= _decomposition_tol else 0
if "symmetric" in mesh:
# Mach-Zehnder interferometers
cmds.append(
Command(
MZgate(np.mod(theta, 2 * np.pi), np.mod(phi, 2 * np.pi)),
(reg[n], reg[m]),
)
)
else:
# Clements style beamsplitters
if not (drop_identity and phi == 0):
cmds.append(Command(Rgate(phi), reg[n]))
if not (drop_identity and theta == 0):
cmds.append(Command(BSgate(theta, 0), (reg[n], reg[m])))
for n, expphi in enumerate(R):
# local phase shifts
q = np.log(expphi).imag if np.abs(expphi - 1) >= _decomposition_tol else 0
if not (drop_identity and q == 0):
cmds.append(Command(Rgate(np.mod(q, 2 * np.pi)), reg[n]))
if BS2 is not None:
# Clements style beamsplitters
for n, m, theta, phi, _ in reversed(BS2):
theta = theta if np.abs(theta) >= _decomposition_tol else 0
phi = phi if np.abs(phi) >= _decomposition_tol else 0
if not (drop_identity and theta == 0):
cmds.append(Command(BSgate(-theta, 0), (reg[n], reg[m])))
if not (drop_identity and phi == 0):
cmds.append(Command(Rgate(-phi), reg[n]))
return cmds
And run the code from step 1 again. This time the output should match as expected.
Why did I change this?
Mathematically speaking I found out that if we use the triangular scheme and want to rebuild the original matrix we have to do the following:
U = T[N-1]^-1 T[N-2]^-1 T[N-3]^-1…T[0]^-1 D
Applying this on a state means that we need to apply D first and then the inverses of T[i]. In the strawberryfields convention a T(theta, phi) = BS(theta, 0) PS(phi). So the inverse would be PS(-phi)* BS(-theta, 0) i.e. changing the order of the operations and the sign of the angle.
This is what i added within the elif mesh == "triangular" block.
Alternatively, you could also change the sf.decomposition.triangular method. Then you would not need to change as much in the _decompose() method from the Interferometer class.
This can be done as follows:
Change the sf.decompositions.triangular method to:
def triangular(V, tol=1e-11):
r"""Triangular decomposition of a unitary matrix due to Reck et al.
See :cite:`reck1994` for more details and :cite:`clements2016` for details on notation.
Args:
V (array[complex]): unitary matrix of size ``n_size``
tol (float): the tolerance used when checking if the matrix is unitary:
:math:`|VV^\dagger-I| \leq` tol
Returns:
tuple[array]: returns a tuple of the form ``(tlist,np.diag(localV), None)``
where:
* ``tlist``: list containing ``[n,m,theta,phi,n_size]`` of the T unitaries needed
* ``localV``: Diagonal unitary applied at the beginning of circuit
"""
localV = V
(nsize, _) = localV.shape
if not np.allclose(V @ V.conj().T, np.identity(nsize), atol=tol, rtol=0):
raise ValueError("The input matrix is not unitary")
tlist = []
for i in range(nsize - 2, -1, -1):
for j in range(i + 1):
tlist.append(nullT(nsize - j - 1, nsize - i - 2, localV))
localV = T(*tlist[-1]) @ localV
return None , np.diag(localV), tlist
Here we only changed the return command from reversed(tlist), np.diag(localV), None to None , np.diag(localV), tlist
In the _decomposition method we then only need to include an if BS1 is not None as follows:
def _decompose(self, reg, **kwargs):
# pylint: disable=too-many-branches
mesh = kwargs.get("mesh", self.mesh)
tol = kwargs.get("tol", self.tol)
drop_identity = kwargs.get("drop_identity", self.drop_identity)
cmds = []
if mesh == "rectangular_compact":
phases = dec.rectangular_compact(self.p[0], rtol=tol, atol=tol)
cmds = _rectangular_compact_cmds(reg, phases)
elif mesh == "triangular_compact":
phases = dec.triangular_compact(self.p[0], rtol=tol, atol=tol)
cmds = _triangular_compact_cmds(reg, phases)
elif mesh == "sun_compact":
parameters, global_phase = dec.sun_compact(self.p[0], rtol=tol, atol=tol)
cmds = _sun_compact_cmds(reg, parameters, global_phase)
elif not self.identity or not drop_identity:
decomp_fn = getattr(dec, mesh)
BS1, R, BS2 = decomp_fn(self.p[0], tol=tol)
if BS1 is not None:
for n, m, theta, phi, _ in BS1:
theta = theta if np.abs(theta) >= _decomposition_tol else 0
phi = phi if np.abs(phi) >= _decomposition_tol else 0
if "symmetric" in mesh:
# Mach-Zehnder interferometers
cmds.append(
Command(
MZgate(np.mod(theta, 2 * np.pi), np.mod(phi, 2 * np.pi)),
(reg[n], reg[m]),
)
)
else:
# Clements style beamsplitters
if not (drop_identity and phi == 0):
cmds.append(Command(Rgate(phi), reg[n]))
if not (drop_identity and theta == 0):
cmds.append(Command(BSgate(theta, 0), (reg[n], reg[m])))
for n, expphi in enumerate(R):
# local phase shifts
q = np.log(expphi).imag if np.abs(expphi - 1) >= _decomposition_tol else 0
if not (drop_identity and q == 0):
cmds.append(Command(Rgate(np.mod(q, 2 * np.pi)), reg[n]))
if BS2 is not None:
# Clements style beamsplitters
for n, m, theta, phi, _ in reversed(BS2):
theta = theta if np.abs(theta) >= _decomposition_tol else 0
phi = phi if np.abs(phi) >= _decomposition_tol else 0
if not (drop_identity and theta == 0):
cmds.append(Command(BSgate(-theta, 0), (reg[n], reg[m])))
if not (drop_identity and phi == 0):
cmds.append(Command(Rgate(-phi), reg[n]))
return cmds
I would say that the second alternative is a bit more clean. You can test both cases and see which one performs better and more reliable.
I hope that this could clarify what I did exactly and I hope that you are able to reproduce this with this instructions. If you need more clarification you can always let me know and I will be happy to help you!
Best wishes,
Jakob
PS: Keep up the good work in the forum. It is really helpful!
This is the new version.