Bug in compile/decompose functions in gaussian backend

[AroosaIjaz]

Issue description

I called the Monte Carlo event probability function which has some issue with the decomposition functions for gaussian backend. I tried to solve the error by changing the tolerance level in the rectangular function in decompositions.py. This did not help and on printing the 'diffn' in the rectangular function, values >1 were also seen. The error looks like this:

----> 7     fv = np.array(similarity.feature_vector_mc(nx.Graph(mol), events, 2, n_mean=8, samples=1000))
      8     t2 = time.time()
      9     fv = fv*(1/np.linalg.norm(fv))
~/Desktop/xanadu/strawberryfields/strawberryfields/apps/similarity.py in feature_vector_mc(graph, event_photon_numbers, max_count_per_mode, n_mean, samples, loss)
    551     return [
    552         prob_event_mc(graph, photons, max_count_per_mode, n_mean, samples, loss)
--> 553         for photons in event_photon_numbers
    554     ]
~/Desktop/xanadu/strawberryfields/strawberryfields/apps/similarity.py in <listcomp>(.0)
    551     return [
    552         prob_event_mc(graph, photons, max_count_per_mode, n_mean, samples, loss)
--> 553         for photons in event_photon_numbers
    554     ]
~/Desktop/xanadu/strawberryfields/strawberryfields/apps/similarity.py in prob_event_mc(graph, photon_number, max_count_per_mode, n_mean, samples, loss)
    451 
    452     eng = sf.LocalEngine(backend="gaussian")
--> 453     result = eng.run(p)
    454 
    455     prob = 0
~/Desktop/xanadu/strawberryfields/strawberryfields/engine.py in run(self, program, args, compile_options, run_options)
    388         eng_run_options = {key: temp_run_options[key] for key in temp_run_options.keys() & eng_run_keys}
    389 
--> 390         result = super()._run(program, args=args, compile_options=compile_options, **eng_run_options)
    391 
    392         modes = temp_run_options["modes"]
~/Desktop/xanadu/strawberryfields/strawberryfields/engine.py in _run(self, program, args, compile_options, **kwargs)
    275             target = self.backend.circuit_spec
    276             if target is not None and p.target != target:
--> 277                 p = p.compile(target, **compile_options)
    278             p.lock()
    279 
~/Desktop/xanadu/strawberryfields/strawberryfields/program.py in compile(self, target, **kwargs)
    509                 )
    510 
--> 511         seq = db.decompose(self.circuit)
    512 
    513         if kwargs.get('warn_connected', True):
~/Desktop/xanadu/strawberryfields/strawberryfields/circuitspecs/circuit_specs.py in decompose(self, seq)
    264                     temp = cmd.op.decompose(cmd.reg, **kwargs)
    265                     # now compile the decomposition
--> 266                     temp = self.decompose(temp)
    267                     compiled.extend(temp)
    268                 except NotImplementedError as err:
~/Desktop/xanadu/strawberryfields/strawberryfields/circuitspecs/circuit_specs.py in decompose(self, seq)
    262                 try:
    263                     kwargs = self.decompositions[op_name]
--> 264                     temp = cmd.op.decompose(cmd.reg, **kwargs)
    265                     # now compile the decomposition
    266                     temp = self.decompose(temp)
~/Desktop/xanadu/strawberryfields/strawberryfields/ops.py in decompose(self, reg, **kwargs)
    418             list[Command]: decomposition as a list of operations acting on specific subsystems
    419         """
--> 420         return self._decompose(reg, **kwargs)
    421 
    422     def _decompose(self, reg, **kwargs):
~/Desktop/xanadu/strawberryfields/strawberryfields/ops.py in _decompose(self, reg, **kwargs)
   1717         if not self.identity or not drop_identity:
   1718             decomp_fn = getattr(dec, mesh)
-> 1719             BS1, R, BS2 = decomp_fn(self.p[0], tol=tol)
   1720 
   1721             for n, m, theta, phi, _ in BS1:
~/Desktop/xanadu/strawberryfields/strawberryfields/decompositions.py in rectangular(V, tol)
    322     diffn = np.linalg.norm(V @ V.conj().T - np.identity(nsize))
    323     if diffn >= tol:
--> 324         raise ValueError("The input matrix is not unitary")
    325 
    326     tilist = []
ValueError: The input matrix is not unitary

• Expected behavior: monte carlo feature vectors outputted

• Actual behavior: ValueError seen

• Reproduces how often: The error might or might not occur when I change the n_mean or the adj matrix.

• System information:

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

Python version: 3.5.6 Platform info: Darwin-18.7.0-x86_64-i386-64bit Installation path: /Users/aroosa/Desktop/xanadu/strawberryfields/strawberryfields Strawberry Fields version: 0.13.0-dev Numpy version: 1.17.4 Scipy version: 1.3.3 SymPy version: 1.5.1 NetworkX version: 2.4 The Walrus version: 0.11.0 Blackbird version: 0.2.3 /opt/anaconda3/envs/graph/lib/python3.5/site-packages/tensorflow/python/framework/dtypes.py:458: FutureWarning:

Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.

/opt/anaconda3/envs/graph/lib/python3.5/site-packages/tensorflow/python/framework/dtypes.py:459: FutureWarning:

Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.

/opt/anaconda3/envs/graph/lib/python3.5/site-packages/tensorflow/python/framework/dtypes.py:460: FutureWarning:

Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.

/opt/anaconda3/envs/graph/lib/python3.5/site-packages/tensorflow/python/framework/dtypes.py:461: FutureWarning:

Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.

/opt/anaconda3/envs/graph/lib/python3.5/site-packages/tensorflow/python/framework/dtypes.py:462: FutureWarning:

Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.

/opt/anaconda3/envs/graph/lib/python3.5/site-packages/tensorflow/python/framework/dtypes.py:465: FutureWarning:

Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.

TensorFlow version: 1.3.0

#### Source code and tracebacks

My code snippet:

```events = [2, 4, 6, 8, 10, 12]```
```fv = np.array(similarity.feature_vector_mc(nx.Graph(test), events, 2, n_mean=8, samples=1000))```

where `test` is the adj matrix.

#### Additional information

I noticed that different functions have different tolerance levels.
Also, I don't think this error comes from Walrus. 

[co9olguy]

Thanks @AroosaIjaz for reporting.

[josh146]

@trbromley, it looks like the apps layer might be generating a non-unitary interferometer?

[trbromley]

With help from @nquesada, we tracked this problem down to the decompositions.takagi() function. This function can return a non-unitary matrix when input matrices that are highly degenerate.

One common case for a degenerate spectrum is when corresponding to a graph, perhaps why @AroosaIjaz noticed this through the graph similarity application in the apps layer.

It's quite an unpredictable bug, i.e., not all degenerate matrices result in non unitary output. However, we managed to find that the adjacency matrix of nx.balanced_tree(2, 4) reliably results in takagi() returning a non-unitary matrix. To accommodate for this, we have updated the takagi function to treat real-valued matrices differently, using a simpler approach that correctly returns a unitary matrix.

I think there may still be a chance that highly degenerate complex-valued matrices can cause a problem, but @nquesada is working on updating takagi() even more.

[co9olguy]

Nice detective work @trbromley @nquesada :shipit:

[AroosaIjaz]

Thank you, guys!