Shannon Decomposition error when synthesizing the GHZ unitary matrix

[ACE07-Sev]

Description of the issue

Greetings there,

Hope all are well. I was trying the quantum_shannon_decomposition() function, and noticed that it doesn't work for the GHZ unitary matrix.

How to reproduce the issue

import cirq
from cirq.transformers.analytical_decompositions import quantum_shannon_decomposition
from scipy.stats import unitary_group

n_qubits = 2
ghz = cirq.Circuit()
qubits = [cirq.NamedQubit(f'q{i}') for i in range(n_qubits)]
ghz.append(cirq.H(qubits[0]))
ghz.append(cirq.CNOT(qubits[0], qubits[1]))

U = ghz.unitary()

qubits = [cirq.NamedQubit(f'q{i}') for i in range(n_qubits)]
circuit = cirq.Circuit(quantum_shannon_decomposition(qubits, U))

assert cirq.approx_eq(U, circuit.unitary(), atol=1e-9)

/home/ace07/Documents/GitHub/QICKIT/.venv/lib/python3.11/site-packages/cirq/transformers/analytical_decompositions/quantum_shannon_decomposition.py:121: RuntimeWarning: divide by zero encountered in scalar divide
  yield ops.ZPowGate(exponent=phi_2 / np.pi, global_shift=phase / phi_2).on(qubit)
/home/ace07/Documents/GitHub/QICKIT/.venv/lib/python3.11/site-packages/cirq/transformers/analytical_decompositions/quantum_shannon_decomposition.py:121: RuntimeWarning: invalid value encountered in scalar divide
  yield ops.ZPowGate(exponent=phi_2 / np.pi, global_shift=phase / phi_2).on(qubit)
/home/ace07/Documents/GitHub/QICKIT/.venv/lib/python3.11/site-packages/cirq/ops/common_gates.py:645: RuntimeWarning: invalid value encountered in scalar multiply
  p = 1j ** (2 * self._exponent * self._global_shift)
/home/ace07/Documents/GitHub/QICKIT/.venv/lib/python3.11/site-packages/cirq/ops/common_gates.py:645: RuntimeWarning: invalid value encountered in scalar multiply
  p = 1j ** (2 * self._exponent * self._global_shift)
  ---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Cell In[7], [line 16](vscode-notebook-cell:?execution_count=7&line=16)
     [13](vscode-notebook-cell:?execution_count=7&line=13) qubits = [cirq.NamedQubit(f'q{i}') for i in range(n_qubits)]
     [14](vscode-notebook-cell:?execution_count=7&line=14) circuit = cirq.Circuit(quantum_shannon_decomposition(qubits, U))
---> [16](vscode-notebook-cell:?execution_count=7&line=16) assert cirq.approx_eq(U, circuit.unitary(), atol=1e-9)

AssertionError:

Cirq version
cirq == 1.4.1

[ACE07-Sev]

If I print the unitary for the circuit, it's actually

array([[nan+nanj, nan+nanj, nan+nanj, nan+nanj],
       [nan+nanj, nan+nanj, nan+nanj, nan+nanj],
       [nan+nanj, nan+nanj, nan+nanj, nan+nanj],
       [nan+nanj, nan+nanj, nan+nanj, nan+nanj]])

Which explains the error.

[ACE07-Sev]

@uzzzzzzz Can you kindly have a look please? I saw that you made the original commit for quantum_shannon_decomposition.py. I have tried a few random circuits, all fail, and give either the AssertionError as the result has nan values only, or they raise the ValueError for non-unitary matrix as .msb_demuxer() creates a non-unitary matrix and sends that to .quantum_shannon_decomposition().

You can try both cases for reproducibility:

• AssertionError: GHZ state, or H(0) Y(1) state with two qubits.
• ValueError: Uniform superposition state, basically Hadamard layer on all qubits.

I understand why the tests have passed before, as the code runs with unitarygroup.rvs(2**num_qubits).

[ACE07-Sev]

So, I implemented this natively on my side, and was able to get past the nan error, but still have the AssertionError problem. It's kind of confusing, it works for certain unitaries, and then for certain unitaries it doesn't.

I tried the following circuit, which is one of the cases that doesn't return the correct value (whilst for randomly generated unitaries, it works fine). A 4 qubit circuit with CX(0, 1) and CX(0, 2).

[NoureldinYosri]

related to #6666