add Quantum Shannon Decomposition (#7907)

* pass tests, lint

* black

* Update qiskit/quantum_info/synthesis/qsd.py

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

* add some documentation

* modify ucry with cz gates

* add release notes. remove unimplmented keyword option.

* remove opt_a2

* remove uneeded import

* add reference to paper

* make demultiplex private.

* remove QuantumShannonDecomposer class.

* Update qiskit/quantum_info/synthesis/qsd.py

Co-authored-by: Ali Javadi-Abhari <[email protected]>

* Update qiskit/quantum_info/synthesis/qsd.py

Co-authored-by: Ali Javadi-Abhari <[email protected]>

Co-authored-by: Matthew Treinish <[email protected]>
Co-authored-by: Ali Javadi-Abhari <[email protected]>
Co-authored-by: mergify[bot] <37929162+mergify[bot]@users.noreply.github.com>

qiskit/extensions/unitary.py

@@ -28,6 +28,7 @@
 from qiskit.quantum_info.operators.predicates import matrix_equal
 from qiskit.quantum_info.operators.predicates import is_unitary_matrix
 from qiskit.quantum_info.synthesis.one_qubit_decompose import OneQubitEulerDecomposer
+from qiskit.quantum_info.synthesis.qsd import qs_decomposition
 from qiskit.quantum_info.synthesis.two_qubit_decompose import two_qubit_cnot_decompose
 from qiskit.extensions.exceptions import ExtensionError
 
@@ -135,10 +136,7 @@ def _define(self):
         elif self.num_qubits == 2:
             self.definition = two_qubit_cnot_decompose(self.to_matrix())
         else:
-            q = QuantumRegister(self.num_qubits, "q")
-            qc = QuantumCircuit(q, name=self.name)
-            qc.append(isometry.Isometry(self.to_matrix(), 0, 0), qargs=q[:])
-            self.definition = qc
+            self.definition = qs_decomposition(self.to_matrix())
 
     def control(self, num_ctrl_qubits=1, label=None, ctrl_state=None):
         """Return controlled version of gate

qiskit/quantum_info/synthesis/qsd.py

@@ -0,0 +1,191 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+"""
+Quantum Shannon Decomposition.
+
+Method is described in arXiv:quant-ph/0406176.
+"""
+import scipy
+import numpy as np
+from qiskit.circuit import QuantumCircuit, QuantumRegister
+from qiskit.quantum_info.synthesis import two_qubit_decompose, one_qubit_decompose
+from qiskit.circuit.library.standard_gates import CXGate
+from qiskit.quantum_info.operators.predicates import is_hermitian_matrix
+from qiskit.extensions.quantum_initializer.uc_pauli_rot import UCPauliRotGate, _EPS
+
+
+def qs_decomposition(mat, opt_a1=True, decomposer_1q=None, decomposer_2q=None):
+    """
+    Decomposes unitary matrix into one and two qubit gates using Quantum Shannon Decomposition.
+
+       ┌───┐               ┌───┐     ┌───┐     ┌───┐
+      ─┤   ├─       ───────┤ Rz├─────┤ Ry├─────┤ Rz├─────
+       │   │    ≃     ┌───┐└─┬─┘┌───┐└─┬─┘┌───┐└─┬─┘┌───┐
+     /─┤   ├─       /─┤   ├──□──┤   ├──□──┤   ├──□──┤   ├
+       └───┘          └───┘     └───┘     └───┘     └───┘
+
+    The number of CX gates generated with the decomposition without optimizations is,
+
+    .. math::
+
+        \frac{9}{16} 4^n - frac{3}{2} 2^n
+
+    If opt_a1=True, the CX count is further reduced by,
+
+    .. math::
+
+        \frac{1}{3} 4^{n - 2} - 1
+
+    This decomposition is described in arXiv:quant-ph/0406176.
+
+    Arguments:
+       mat (ndarray): unitary matrix to decompose
+       opt_a1 (bool): whether to try optimization A.1 from Shende. This should eliminate 1 cnot
+          per call. If True CZ gates are left in the output. If desired these can be further decomposed
+          to CX.
+       decomposer_1q (None or Object): optional 1Q decomposer. If None, uses
+          :class:`~qiskit.quantum_info.synthesis.one_qubit_decomposer.OneQubitEulerDecomser`
+       decomposer_2q (None or Object): optional 2Q decomposer. If None, uses
+          :class:`~qiskit.quantum_info.synthesis.two_qubit_decomposer.TwoQubitBasisDecomposer`
+          with CXGate.
+
+    Return:
+       QuantumCircuit: Decomposed quantum circuit.
+    """
+    dim = mat.shape[0]
+    nqubits = int(np.log2(dim))
+    if np.allclose(np.identity(dim), mat):
+        return QuantumCircuit(nqubits)
+    if dim == 2:
+        if decomposer_1q is None:
+            decomposer_1q = one_qubit_decompose.OneQubitEulerDecomposer()
+        circ = decomposer_1q(mat)
+    elif dim == 4:
+        if decomposer_2q is None:
+            decomposer_2q = two_qubit_decompose.TwoQubitBasisDecomposer(CXGate())
+        circ = decomposer_2q(mat)
+    else:
+        qr = QuantumRegister(nqubits)
+        circ = QuantumCircuit(qr)
+        dim_o2 = dim // 2
+        # perform cosine-sine decomposition
+        (u1, u2), vtheta, (v1h, v2h) = scipy.linalg.cossin(mat, separate=True, p=dim_o2, q=dim_o2)
+        # left circ
+        left_circ = _demultiplex(v1h, v2h, opt_a1=opt_a1)
+        circ.append(left_circ.to_instruction(), qr)
+        # middle circ
+        if opt_a1:
+            nangles = len(vtheta)
+            half_size = nangles // 2
+            # get UCG in terms of CZ
+            circ_cz = _get_ucry_cz(nqubits, (2 * vtheta).tolist())
+            circ.append(circ_cz.to_instruction(), range(nqubits))
+            # merge final cz with right-side generic multiplexer
+            u2[:, half_size:] = np.negative(u2[:, half_size:])
+        else:
+            circ.ucry((2 * vtheta).tolist(), qr[:-1], qr[-1])
+        # right circ
+        right_circ = _demultiplex(u1, u2, opt_a1=opt_a1)
+        circ.append(right_circ.to_instruction(), qr)
+
+    return circ
+
+
+def _demultiplex(um0, um1, opt_a1=False):
+    """decomposes a generic multiplexer.
+
+          ────□────
+           ┌──┴──┐
+         /─┤     ├─
+           └─────┘
+
+    represented by the block diagonal matrix
+
+            ┏         ┓
+            ┃ um0     ┃
+            ┃     um1 ┃
+            ┗         ┛
+
+    to
+               ┌───┐
+        ───────┤ Rz├──────
+          ┌───┐└─┬─┘┌───┐
+        /─┤ w ├──□──┤ v ├─
+          └───┘     └───┘
+
+    where v and w are general unitaries determined from decomposition.
+
+    Args:
+       um0 (ndarray): applied if MSB is 0
+       um1 (ndarray): applied if MSB is 1
+       opt_a1 (bool): whether to try optimization A.1 from Shende. This should elliminate 1 cnot
+          per call. If True CZ gates are left in the output. If desired these can be further decomposed
+
+    Returns:
+        QuantumCircuit: decomposed circuit
+    """
+    dim = um0.shape[0] + um1.shape[0]  # these should be same dimension
+    nqubits = int(np.log2(dim))
+    um0um1 = um0 @ um1.T.conjugate()
+    if is_hermitian_matrix(um0um1):
+        eigvals, vmat = scipy.linalg.eigh(um0um1)
+    else:
+        evals, vmat = scipy.linalg.schur(um0um1, output="complex")
+        eigvals = evals.diagonal()
+    dvals = np.lib.scimath.sqrt(eigvals)
+    dmat = np.diag(dvals)
+    wmat = dmat @ vmat.T.conjugate() @ um1
+
+    circ = QuantumCircuit(nqubits)
+
+    # left gate
+    left_gate = qs_decomposition(wmat, opt_a1=opt_a1).to_instruction()
+    circ.append(left_gate, range(nqubits - 1))
+
+    # multiplexed Rz
+    angles = 2 * np.angle(np.conj(dvals))
+    circ.ucrz(angles.tolist(), list(range(nqubits - 1)), [nqubits - 1])
+
+    # right gate
+    right_gate = qs_decomposition(vmat, opt_a1=opt_a1).to_instruction()
+    circ.append(right_gate, range(nqubits - 1))
+
+    return circ
+
+
+def _get_ucry_cz(nqubits, angles):
+    """
+    Get uniformally controlled Ry gate in in CZ-Ry as in UCPauliRotGate.
+    """
+    nangles = len(angles)
+    qc = QuantumCircuit(nqubits)
+    q_controls = qc.qubits[:-1]
+    q_target = qc.qubits[-1]
+    if not q_controls:
+        if np.abs(angles[0]) > _EPS:
+            qc.ry(angles[0], q_target)
+    else:
+        angles = angles.copy()
+        UCPauliRotGate._dec_uc_rotations(angles, 0, len(angles), False)
+        for (i, angle) in enumerate(angles):
+            if np.abs(angle) > _EPS:
+                qc.ry(angle, q_target)
+            if not i == len(angles) - 1:
+                binary_rep = np.binary_repr(i + 1)
+                q_contr_index = len(binary_rep) - len(binary_rep.rstrip("0"))
+            else:
+                # Handle special case:
+                q_contr_index = len(q_controls) - 1
+            # leave off last CZ for merging with adjacent UCG
+            if i < nangles - 1:
+                qc.cz(q_controls[q_contr_index], q_target)
+    return qc

qiskit/transpiler/passes/synthesis/unitary_synthesis.py

@@ -22,10 +22,10 @@
 from qiskit.transpiler.basepasses import TransformationPass
 from qiskit.transpiler.exceptions import TranspilerError
 from qiskit.dagcircuit.dagcircuit import DAGCircuit
-from qiskit.extensions.quantum_initializer import isometry
 from qiskit.quantum_info.synthesis import one_qubit_decompose
 from qiskit.quantum_info.synthesis.xx_decompose import XXDecomposer
 from qiskit.quantum_info.synthesis.two_qubit_decompose import TwoQubitBasisDecomposer
+from qiskit.quantum_info.synthesis.qsd import qs_decomposition
 from qiskit.circuit.parameter import Parameter
 from qiskit.circuit.library.standard_gates import (
     iSwapGate,
@@ -539,7 +539,7 @@ def run(self, unitary, **options):
                 preferred_direction,
             )
         else:
-            synth_dag = circuit_to_dag(isometry.Isometry(unitary, 0, 0).definition)
+            synth_dag = circuit_to_dag(qs_decomposition(unitary))
 
         return synth_dag, wires
 

releasenotes/notes/quantum_shannon_decomp-facaa362a3ca8ba3.yaml

@@ -0,0 +1,13 @@
+---
+features:
+  - |
+    Added a new module to apply Quantum Shannon Decomposition of
+    arbitrary unitaries. This functionality replaces the previous
+    usage of the isometry module in the unitary synthesis transpiler
+    pass as well as when adding unitaries to a circuit using a
+    UnitaryGate.
+upgrade:
+  - |
+    The Quantum Shannon Decomposition uses about half the cnot gates as the
+    isometry implementation when decomposing unitary matrices of greater than
+    two qubits.

test/python/quantum_info/test_synthesis.py

@@ -17,10 +17,12 @@
 import logging
 from test import combine
 
-from ddt import ddt
 import numpy as np
+import scipy
+import scipy.stats
+from ddt import ddt, data
 
-from qiskit import execute, QiskitError
+from qiskit import execute, QiskitError, transpile
 from qiskit.circuit import QuantumCircuit, QuantumRegister
 from qiskit.converters import dag_to_circuit, circuit_to_dag
 from qiskit.extensions import UnitaryGate
@@ -72,6 +74,7 @@
 )
 
 from qiskit.quantum_info.synthesis.ion_decompose import cnot_rxx_decompose
+import qiskit.quantum_info.synthesis.qsd as qsd
 from qiskit.test import QiskitTestCase
 
 
@@ -349,14 +352,14 @@ def check_oneq_special_cases(
         """Check OneQubitEulerDecomposer produces the expected gates"""
         decomposer = OneQubitEulerDecomposer(basis)
         circ = decomposer(target, simplify=True)
-        data = Operator(circ).data
-        maxdist = np.max(np.abs(target.data - data))
-        trace = np.trace(data.T.conj() @ target)
+        cdata = Operator(circ).data
+        maxdist = np.max(np.abs(target.data - cdata))
+        trace = np.trace(cdata.T.conj() @ target)
         self.assertLess(
             np.abs(maxdist),
             tolerance,
             f"Worst case distance: {maxdist}, trace: {trace}\n"
-            f"Target:\n{target}\nActual:\n{data}\n{circ}",
+            f"Target:\n{target}\nActual:\n{cdata}\n{circ}",
         )
         if expected_gates is not None:
             self.assertDictEqual(dict(circ.count_ops()), expected_gates, f"Circuit:\n{circ}")
@@ -1360,5 +1363,103 @@ def test_decompose_two_qubit_product_gate_not_product(self):
         self.assertIn("decomposition failed", exc.exception.message)
 
 
+@ddt
+class TestQuantumShannonDecomposer(QiskitTestCase):
+    """
+    Test Quantum Shannon Decomposition.
+    """
+
+    def setUp(self):
+        super().setUp()
+        seed = (hash(self.id())) % 10000
+        np.random.seed(seed)
+        self.qsd = qsd.qs_decomposition
+
+    def _get_lower_cx_bound(self, n):
+        return 1 / 4 * (4**n - 3 * n - 1)
+
+    def _qsd_l2_cx_count(self, n):
+        """expected unoptimized cnot count for down to 2q"""
+        return 9 / 16 * 4**n - 3 / 2 * 2**n
+
+    def _qsd_l2_a1_mod(self, n):
+        return (4 ** (n - 2) - 1) // 3
+
+    def _qsd_l2_a2_mod(self, n):
+        return 4 ** (n - 1) - 1
+
+    @data(*list(range(1, 5)))
+    def test_random_decomposition_l2_no_opt(self, nqubits):
+        """test decomposition of random SU(n) down to 2 qubits without optimizations."""
+        dim = 2**nqubits
+        mat = scipy.stats.unitary_group.rvs(dim)
+        circ = self.qsd(mat, opt_a1=False)
+        ccirc = transpile(circ, basis_gates=["u", "cx"], optimization_level=0)
+        self.assertTrue(np.allclose(mat, Operator(ccirc).data))
+        if nqubits > 1:
+            self.assertLessEqual(ccirc.count_ops().get("cx"), self._qsd_l2_cx_count(nqubits))
+        else:
+            self.assertEqual(sum(ccirc.count_ops().values()), 1)
+
+    @data(*list(range(1, 5)))
+    def test_random_decomposition_l2_a1_opt(self, nqubits):
+        """test decomposition of random SU(n) down to 2 qubits with 'a1' optimization."""
+        dim = 2**nqubits
+        mat = scipy.stats.unitary_group.rvs(dim)
+        circ = self.qsd(mat, opt_a1=True)
+        ccirc = transpile(circ, basis_gates=["u", "cx"], optimization_level=0)
+        self.assertTrue(np.allclose(mat, Operator(ccirc).data))
+        if nqubits > 1:
+            expected_cx = self._qsd_l2_cx_count(nqubits) - self._qsd_l2_a1_mod(nqubits)
+            self.assertLessEqual(ccirc.count_ops().get("cx"), expected_cx)
+
+    def test_SO3_decomposition_l2_a1_opt(self):
+        """test decomposition of random So(3) down to 2 qubits with 'a1' optimization."""
+        nqubits = 3
+        dim = 2**nqubits
+        mat = scipy.stats.ortho_group.rvs(dim)
+        circ = self.qsd(mat, opt_a1=True)
+        ccirc = transpile(circ, basis_gates=["u", "cx"], optimization_level=0)
+        self.assertTrue(np.allclose(mat, Operator(ccirc).data))
+        expected_cx = self._qsd_l2_cx_count(nqubits) - self._qsd_l2_a1_mod(nqubits)
+        self.assertLessEqual(ccirc.count_ops().get("cx"), expected_cx)
+
+    def test_identity_decomposition(self):
+        """Test decomposition on identity matrix"""
+        nqubits = 3
+        dim = 2**nqubits
+        mat = np.identity(dim)
+        circ = self.qsd(mat, opt_a1=True)
+        self.assertTrue(np.allclose(mat, Operator(circ).data))
+        self.assertEqual(sum(circ.count_ops().values()), 0)
+
+    @data(*list(range(1, 4)))
+    def test_diagonal(self, nqubits):
+        """Test decomposition on diagonal -- qsd is not optimal"""
+        dim = 2**nqubits
+        mat = np.diag(np.exp(1j * np.random.normal(size=dim)))
+        circ = self.qsd(mat, opt_a1=True)
+        ccirc = transpile(circ, basis_gates=["u", "cx"], optimization_level=0)
+        self.assertTrue(np.allclose(mat, Operator(ccirc).data))
+        if nqubits > 1:
+            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(2, 4)))
+    def test_hermitian(self, nqubits):
+        """Test decomposition on hermitian -- qsd is not optimal"""
+        # better might be (arXiv:1405.6741)
+        dim = 2**nqubits
+        umat = scipy.stats.unitary_group.rvs(dim)
+        dmat = np.diag(np.exp(1j * np.random.normal(size=dim)))
+        mat = umat.T.conjugate() @ dmat @ umat
+        circ = self.qsd(mat, opt_a1=True)
+        ccirc = transpile(circ, basis_gates=["u", "cx"], optimization_level=0)
+        self.assertTrue(np.allclose(mat, Operator(ccirc).data))
+        if nqubits > 1:
+            expected_cx = self._qsd_l2_cx_count(nqubits) - self._qsd_l2_a1_mod(nqubits)
+            self.assertLessEqual(ccirc.count_ops().get("cx"), expected_cx)
+
+
 if __name__ == "__main__":
     unittest.main()