add test of block extration wrt qubit

qiskit/synthesis/unitary/qsd.py

@@ -19,7 +19,7 @@
 from typing import Callable
 import scipy
 import numpy as np
-from qiskit.circuit.quantumcircuit import QuantumCircuit, QuantumRegister
+from qiskit.circuit.quantumcircuit import QuantumCircuit, QuantumRegister, Qubit
 from qiskit.synthesis.two_qubit import (
     TwoQubitBasisDecomposer,
     two_qubit_decompose,
@@ -121,6 +121,13 @@ def decomp_2q(mat):
                 decomposer_2q = TwoQubitBasisDecomposer(CXGate())
         circ = decomposer_2q(mat)
     else:
+        # check whether matrix is equivalent to block diagonal wrt ctrl_index
+        if opt_a2 is False:
+            for ctrl_index in range(nqubits):
+                um00, um11, um01, um10 = _extract_multiplex_blocks(mat, ctrl_index)
+                if _off_diagonals_are_zero(um01, um10):
+                    return _demultiplex(um01, um11, opt_a1=opt_a1, opt_a2=opt_a2,
+                                        _depth=_depth, _ctrl_index=ctrl_index)
         qr = QuantumRegister(nqubits)
         circ = QuantumCircuit(qr)
         dim_o2 = dim // 2
@@ -150,7 +157,7 @@ def decomp_2q(mat):
     return circ
 
 
-def _demultiplex(um0, um1, opt_a1=False, opt_a2=False, *, _depth=0):
+def _demultiplex(um0, um1, opt_a1=False, opt_a2=False, *, _depth=0, _ctrl_index=0):
     """Decompose a generic multiplexer.
 
           ────□────
@@ -183,6 +190,7 @@ def _demultiplex(um0, um1, opt_a1=False, opt_a2=False, *, _depth=0):
           unitaries into a diagonal gate and a two cx unitary and reduces overall cx count by
           4^(n-2) - 1.
        _depth (int): This is an internal variable to track the recursion depth.
+       _ctrl_index (int): The index wrt which um0 and um1 are controlled.
 
     Returns:
         QuantumCircuit: decomposed circuit
@@ -199,24 +207,27 @@ def _demultiplex(um0, um1, opt_a1=False, opt_a2=False, *, _depth=0):
     dmat = np.diag(dvals)
     wmat = dmat @ vmat.T.conjugate() @ um1
 
-    circ = QuantumCircuit(nqubits)
+    qr = [Qubit() for _ in range(nqubits)]
+    circ = QuantumCircuit(qr)
 
     # left gate
     left_gate = qs_decomposition(
         wmat, opt_a1=opt_a1, opt_a2=opt_a2, _depth=_depth + 1
     ).to_instruction()
-    circ.append(left_gate, range(nqubits - 1))
+
+    qr_sub = [qr[i] for i in range(nqubits) if i != _ctrl_index]
+    circ.append(left_gate, qr_sub)
 
     # multiplexed Rz
     angles = 2 * np.angle(np.conj(dvals))
     ucrz = UCRZGate(angles.tolist())
-    circ.append(ucrz, [nqubits - 1] + list(range(nqubits - 1)))
+    circ.append(ucrz, [qr[_ctrl_index]] + qr_sub)
 
     # right gate
     right_gate = qs_decomposition(
         vmat, opt_a1=opt_a1, opt_a2=opt_a2, _depth=_depth + 1
     ).to_instruction()
-    circ.append(right_gate, range(nqubits - 1))
+    circ.append(right_gate, qr_sub)
 
     return circ
 
@@ -286,3 +297,42 @@ def _apply_a2(circ):
     qc3 = two_qubit_decompose.two_qubit_cnot_decompose(mat2)
     ccirc.data[ind2] = ccirc.data[ind2].replace(operation=qc3.to_gate())
     return ccirc
+
+def _extract_multiplex_blocks(U, k):
+    """
+    A block diagonal gate is represented as 
+    """
+    dim = U.shape[0]
+    N = dim.bit_length() - 1
+    halfdim = dim // 2
+
+    U_tensor = U.reshape((2,) * N + (2,) * N)
+
+    # Move qubit k to top
+    if k != 0:
+        U_tensor = np.moveaxis(U_tensor, k,   0)
+        U_tensor = np.moveaxis(U_tensor, k+N, N)
+
+    # reshape for extraction
+    U_4d = U_tensor.reshape(2, halfdim, 2, halfdim)
+    # block for qubit k = |0>
+    um00 = U_4d[0, :, 0, :]
+    # block for qubit k = |1>
+    um11 = U_4d[1, :, 1, :]
+    # off diagonal blocks
+    um01 = U_4d[0, :, 1, :]
+    um10 = U_4d[1, :, 0, :]
+    return um00, um11, um01, um10
+
+def _off_diagonals_are_zero(um01, um10, atol=1e-12):
+    """
+    Checks whether off-diagonal blocks are zero.
+    Args:
+       um01 (ndarray): upper right block
+       um10 (ndarray): lower left block
+       atol (float): absolute tolerance
+    Returns:
+       bool: whether both blocks are zero within tolerance
+    """
+    return (np.allclose(um01, 0, atol=atol) and
+            np.allclose(um10, 0, atol=atol))

test/python/synthesis/test_synthesis.py

@@ -17,6 +17,7 @@
 import contextlib
 import logging
 import math
+import itertools
 import numpy as np
 import scipy
 import scipy.stats
@@ -54,6 +55,7 @@
     RYGate,
     RZGate,
     UnitaryGate,
+    UCGate,
 )
 from qiskit.quantum_info.operators import Operator
 from qiskit.quantum_info.random import random_unitary
@@ -1769,6 +1771,87 @@ def test_a2_opt_single_2q(self):
         except UnboundLocalError as uerr:
             self.fail(str(uerr))
 
+    def _create_random_multiplexed_gate(self, num_qubits):
+        num_blocks = 2
+        blocks = [scipy.stats.unitary_group.rvs(2**(num_qubits-1)) for _ in range(num_blocks)]
+        mat = scipy.linalg.block_diag(*blocks)  # control on "top" qubit
+        multiplexed_gate = UnitaryGate(mat)
+        return multiplexed_gate, mat
+
+    def test_tensor_block_uc_2q(self):
+        np.set_printoptions(linewidth=250, precision=2, suppress=True)
+        num_qubits = 2
+        gate, mat = self._create_random_multiplexed_gate(num_qubits)
+        print(f'original matrix:\n', mat)
+        for layout in itertools.permutations(range(num_qubits)):
+            print('='*30)
+            # create gate with "control" on different qubits
+            qc = QuantumCircuit(num_qubits)
+            qc.append(gate, layout)
+            hidden_mat = Operator(qc).data
+            print(f'layout = {layout}')
+            print(hidden_mat)
+            for j in range(num_qubits):
+                um00, um11, um01, um10 = qsd._extract_multiplex_blocks(
+                    hidden_mat, k=j)
+                # Check the off-diagonal
+                is_mult = qsd._off_diagonals_are_zero(um01, um10)
+                print(f"Is multiplexed wrt qubit {j}? {is_mult}")
+                if is_mult:
+                    qc_uc = QuantumCircuit(num_qubits)
+                    ucgate = UCGate([um00, um11])
+                    ucgate.definition
+                    qc_uc.append(ucgate, layout)
+                    uc_op = Operator(qc_uc)
+                    self.assertEqual(Operator(qc), uc_op)
+
+    def _get_multiplex_matrix(self, um00, um11, k):
+        """form matrix multiplexed wrt qubit k"""
+        halfdim = um00.shape[0]
+        dim = 2 * halfdim
+        N = halfdim.bit_length()
+        Ure4d = np.zeros((2, halfdim, 2, halfdim), dtype=complex)
+        Ure4d[0,:,0,:] = um00
+        Ure4d[1,:,1,:] = um11
+        UreNd = Ure4d.reshape((2,)*N + (2,)*N)
+        UreNd = np.moveaxis(UreNd, N, k+N)
+        UreNd = np.moveaxis(UreNd, 0, k)
+        Ure = UreNd.reshape(dim, dim)
+        return Ure
+
+    def test_tensor_block_uc_3q(self):
+        """Create 3q gate with multiplexed controls"""
+        np.set_printoptions(linewidth=250, precision=2, suppress=True)
+        num_qubits = 3
+        gate, mat = self._create_random_multiplexed_gate(num_qubits)
+        print(f'original matrix:\n', mat)
+        for layout in itertools.permutations(range(num_qubits)):
+            print('='*30)
+            # create gate with "control" on different qubits
+            qc = QuantumCircuit(num_qubits)
+            qc.append(gate, layout)
+            hidden_mat = Operator(qc).data
+            print(f'layout = {layout}')
+            print(hidden_mat)
+            for j in range(num_qubits):
+                um00, um11, um01, um10 = qsd._extract_multiplex_blocks(hidden_mat, k=j)
+                # Check the off-diagonal
+                is_mult = qsd._off_diagonals_are_zero(um01, um10)
+                print(f"Is multiplexed wrt qubit {j}? {is_mult}")
+                if is_mult:
+                    qc_uc = QuantumCircuit(num_qubits)
+                    uc_mat = self._get_multiplex_matrix(um00, um11, j)
+                    uc_gate = UnitaryGate(uc_mat)
+                    qc_uc.append(uc_gate, range(num_qubits))
+                    uc_op = Operator(qc_uc)
+                    if Operator(qc) == uc_op:
+                        print(f'VALIDATED: {layout}, index={j}')
+                    else:
+                        print(f'NOT VALID: {layout}, index={j}')
+                        print(uc_mat)
+                        print('um01:\n', um01)
+                        print('um10:\n', um10)
+
 
 class TestTwoQubitDecomposeUpToDiagonal(QiskitTestCase):
     """test TwoQubitDecomposeUpToDiagonal class"""