Bloch-Messiah returns incorrect results

[nquesada]

Before posting a bug report

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

Expected behavior

Consider

import numpy as np
from thewalrus.random import random_symplectic
from thewalrus.quantum import is_symplectic
from strawberryfields.decompositions import bloch_messiah

r = np.array([1,1,1,1,2,2,0,0,0,0,0,0])
Sdiag= np.diag(np.concatenate([np.exp(r), np.exp(-r)]))
num_modes = len(r)
O1 = random_symplectic(num_modes, passive=True)
O2 = random_symplectic(num_modes, passive=True)
S = O1 @ Sdiag @ O2
assert is_symplectic(S)

# Do the Bloch-Messiah from SF
from strawberryfields.decompositions import bloch_messiah
O1, S, O2 = bloch_messiah(S)

assert is_symplectic(Sdiag)
assert is_symplectic(O1)
assert is_symplectic(O2)

Actual behavior

The returned matrices O1 and O2 in the example above are not symplectic.

Reproduces how often

Always

System information

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

Python version:            3.9.7
Platform info:             
Installation path:         
Strawberry Fields version: 0.22.0-dev
Numpy version:             1.21.5
Scipy version:             1.8.0
SymPy version:             1.10
NetworkX version:          2.6.3
The Walrus version:        0.20.0-dev
Blackbird version:         0.4.0
XCC version:               0.1.2
TensorFlow version:        2.6.2

Source code

No response

Tracebacks

No response

Additional information

A fixed Bloch-Messiah can be found in https://github.com/polyquantique/thewalrus/pull/2

[josh146]

Thanks for catching this @nquesada 😬

[sduquemesa]

I haven't been able to find the root cause of this but it also happens for r = np.array([1,0,0]) or any other array with more than one zero. Any ideas?