fix 1q matrix bug in Quantum Shannon Decomposer (#10126)

* fix 1q bug

* formatting

* restrict 2q gates from apply_a2

* Add check that decomposition includes qsd2q gates before optimizing them

Co-authored-by: jsmallz333 <[email protected]>

* Apply suggestions from code review

Co-authored-by: Matthew Treinish <[email protected]>

* Add release note for fix to issue 10036

* avoid creating unitary gate if initial gate is 2q

* remove debug code

* Change "certain" to "trivial"

This was requested in a review comment.

---------

Co-authored-by: John Lapeyre <[email protected]>
Co-authored-by: jsmallz333 <[email protected]>
Co-authored-by: Matthew Treinish <[email protected]>
(cherry picked from commit b831bcf)

qiskit/quantum_info/synthesis/qsd.py

@@ -81,7 +81,7 @@ def qs_decomposition(
         circ = decomposer_1q(mat)
     elif dim == 4:
         if decomposer_2q is None:
-            if opt_a2:
+            if opt_a2 and _depth > 0:
                 from qiskit.extensions.unitary import UnitaryGate  # pylint: disable=cyclic-import
 
                 def decomp_2q(mat):
@@ -118,7 +118,7 @@ def decomp_2q(mat):
         right_circ = _demultiplex(u1, u2, opt_a1=opt_a1, opt_a2=opt_a2, _depth=_depth)
         circ.append(right_circ.to_instruction(), qr)
 
-    if opt_a2 and _depth == 0:
+    if opt_a2 and _depth == 0 and dim > 4:
         return _apply_a2(circ)
     return circ
 
@@ -237,6 +237,8 @@ def _apply_a2(circ):
         instr, _, _ = instr_context
         if instr.name == "qsd2q":
             ind2q.append(i)
+    if not ind2q:
+        return ccirc
     # rolling over diagonals
     ind2 = None  # lint
     for ind1, ind2 in zip(ind2q[0:-1:], ind2q[1::]):

releasenotes/notes/fix-1q-matrix-bug-in-quantum-shannon-decomposer-c99ce6509f03715b.yaml

@@ -0,0 +1,7 @@
+---
+fixes:
+  - |
+    When :func:`~qiskit.quantum_info.synthesis.qs_decomposition`, which does quantum Shannon
+    decomposition, was called on trivial numeric unitaries that do not benefit from this
+    decomposition, an unexpected error was raised. With this fix, such unitaries are detected and
+    the equivalent circuit is returned without performing Shannon decomposition.

test/python/quantum_info/test_synthesis.py

@@ -1530,17 +1530,23 @@ def test_hermitian(self, nqubits):
             expected_cx = self._qsd_l2_cx_count(nqubits) - self._qsd_l2_a1_mod(nqubits)
             self.assertLessEqual(ccirc.count_ops().get("cx"), expected_cx)
 
-    @data(*list(range(3, 6)))
+    @data(*list(range(1, 6)))
     def test_opt_a1a2(self, nqubits):
         """Test decomposition with both optimization a1 and a2 from shende2006"""
         dim = 2**nqubits
         umat = scipy.stats.unitary_group.rvs(dim, random_state=1224)
         circ = self.qsd(umat, opt_a1=True, opt_a2=True)
         ccirc = transpile(circ, basis_gates=["u", "cx"], optimization_level=0)
         self.assertTrue(Operator(umat) == Operator(ccirc))
-        self.assertEqual(
-            ccirc.count_ops().get("cx"), (23 / 48) * 4**nqubits - (3 / 2) * 2**nqubits + 4 / 3
-        )
+        if nqubits > 2:
+            self.assertEqual(
+                ccirc.count_ops().get("cx"),
+                (23 / 48) * 4**nqubits - (3 / 2) * 2**nqubits + 4 / 3,
+            )
+        elif nqubits == 1:
+            self.assertEqual(ccirc.count_ops().get("cx", 0), 0)
+        elif nqubits == 2:
+            self.assertLessEqual(ccirc.count_ops().get("cx", 0), 3)
 
 
 class TestTwoQubitDecomposeUpToDiagonal(QiskitTestCase):