Add Clifford.from_matrix (Qiskit#9475)

* Add Clifford.from_matrix

* Add tests

* Add reno

* Faster Clifford.from_matrix O(16^n)->O(4^n)

* Add infinite recursion test case

* Change to try append with definition first

* Add u gate handling for speed

* Improve implematation following review comments

* Add Clifford.from_operator

* Update reno

* Lint

* more accurate reno

---------

Co-authored-by: Ikko Hamamura <[email protected]>

qiskit/quantum_info/operators/symplectic/clifford.py

@@ -12,6 +12,7 @@
 """
 Clifford operator class.
 """
+import functools
 import itertools
 import re
 
@@ -552,10 +553,50 @@ def to_matrix(self):
         """Convert operator to Numpy matrix."""
         return self.to_operator().data
 
+    @classmethod
+    def from_matrix(cls, matrix):
+        """Create a Clifford from a unitary matrix.
+
+        Note that this function takes exponentially long time w.r.t. the number of qubits.
+
+        Args:
+            matrix (np.array): A unitary matrix representing a Clifford to be converted.
+
+        Returns:
+            Clifford: the Clifford object for the unitary matrix.
+
+        Raises:
+            QiskitError: if the input is not a Clifford matrix.
+        """
+        tableau = cls._unitary_matrix_to_tableau(matrix)
+        if tableau is None:
+            raise QiskitError("Non-Clifford matrix is not convertible")
+        return cls(tableau)
+
     def to_operator(self):
         """Convert to an Operator object."""
         return Operator(self.to_instruction())
 
+    @classmethod
+    def from_operator(cls, operator):
+        """Create a Clifford from a operator.
+
+        Note that this function takes exponentially long time w.r.t. the number of qubits.
+
+        Args:
+            operator (Operator): An operator representing a Clifford to be converted.
+
+        Returns:
+            Clifford: the Clifford object for the operator.
+
+        Raises:
+            QiskitError: if the input is not a Clifford operator.
+        """
+        tableau = cls._unitary_matrix_to_tableau(operator.to_matrix())
+        if tableau is None:
+            raise QiskitError("Non-Clifford operator is not convertible")
+        return cls(tableau)
+
     def to_circuit(self):
         """Return a QuantumCircuit implementing the Clifford.
 
@@ -604,9 +645,9 @@ def from_circuit(circuit):
         # Initialize an identity Clifford
         clifford = Clifford(np.eye(2 * circuit.num_qubits), validate=False)
         if isinstance(circuit, QuantumCircuit):
-            _append_circuit(clifford, circuit)
+            clifford = _append_circuit(clifford, circuit)
         else:
-            _append_operation(clifford, circuit)
+            clifford = _append_operation(clifford, circuit)
         return clifford
 
     @staticmethod
@@ -662,7 +703,7 @@ def from_label(label):
         num_qubits = len(label)
         op = Clifford(np.eye(2 * num_qubits, dtype=bool))
         for qubit, char in enumerate(reversed(label)):
-            _append_operation(op, label_gates[char], qargs=[qubit])
+            op = _append_operation(op, label_gates[char], qargs=[qubit])
         return op
 
     def to_labels(self, array=False, mode="B"):
@@ -849,6 +890,95 @@ def _from_label(label):
         symp[-1] = phase
         return symp
 
+    @staticmethod
+    def _pauli_matrix_to_row(mat, num_qubits):
+        """Generate a binary vector (a row of tableau representation) from a Pauli matrix.
+        Return None if the non-Pauli matrix is supplied."""
+        # pylint: disable=too-many-return-statements
+
+        def find_one_index(x, decimals=6):
+            indices = np.where(np.round(np.abs(x), decimals) == 1)
+            return indices[0][0] if len(indices[0]) == 1 else None
+
+        def bitvector(n, num_bits):
+            return np.array([int(digit) for digit in format(n, f"0{num_bits}b")], dtype=bool)[::-1]
+
+        # compute x-bits
+        xint = find_one_index(mat[0, :])
+        if xint is None:
+            return None
+        xbits = bitvector(xint, num_qubits)
+
+        # extract non-zero elements from matrix (rounded to 1, -1, 1j or -1j)
+        entries = np.empty(len(mat), dtype=complex)
+        for i, row in enumerate(mat):
+            index = find_one_index(row)
+            if index is None:
+                return None
+            expected = xint ^ i
+            if index != expected:
+                return None
+            entries[i] = np.round(mat[i, index])
+
+        # compute z-bits
+        zbits = np.empty(num_qubits, dtype=bool)
+        for k in range(num_qubits):
+            sign = np.round(entries[2**k] / entries[0])
+            if sign == 1:
+                zbits[k] = False
+            elif sign == -1:
+                zbits[k] = True
+            else:
+                return None
+
+        # compute phase
+        phase = None
+        num_y = sum(xbits & zbits)
+        positive_phase = (-1j) ** num_y
+        if entries[0] == positive_phase:
+            phase = False
+        elif entries[0] == -1 * positive_phase:
+            phase = True
+        if phase is None:
+            return None
+
+        # validate all non-zero elements
+        coef = ((-1) ** phase) * positive_phase
+        ivec, zvec = np.ones(2), np.array([1, -1])
+        expected = coef * functools.reduce(np.kron, [zvec if z else ivec for z in zbits[::-1]])
+        if not np.allclose(entries, expected):
+            return None
+
+        return np.hstack([xbits, zbits, phase])
+
+    @staticmethod
+    def _unitary_matrix_to_tableau(matrix):
+        # pylint: disable=invalid-name
+        num_qubits = int(np.log2(len(matrix)))
+
+        stab = np.empty((num_qubits, 2 * num_qubits + 1), dtype=bool)
+        for i in range(num_qubits):
+            label = "I" * (num_qubits - i - 1) + "X" + "I" * i
+            Xi = Operator.from_label(label).to_matrix()
+            target = matrix @ Xi @ np.conj(matrix).T
+            row = Clifford._pauli_matrix_to_row(target, num_qubits)
+            if row is None:
+                return None
+            stab[i] = row
+
+        destab = np.empty((num_qubits, 2 * num_qubits + 1), dtype=bool)
+        for i in range(num_qubits):
+            label = "I" * (num_qubits - i - 1) + "Z" + "I" * i
+            Zi = Operator.from_label(label).to_matrix()
+            target = matrix @ Zi @ np.conj(matrix).T
+            row = Clifford._pauli_matrix_to_row(target, num_qubits)
+            if row is None:
+                return None
+            destab[i] = row
+
+        tableau = np.vstack([stab, destab])
+        return tableau
+
 
 # Update docstrings for API docs
 generate_apidocs(Clifford)

qiskit/quantum_info/operators/symplectic/clifford_circuits.py

@@ -12,19 +12,22 @@
 """
 Circuit simulation for the Clifford class.
 """
+import copy
+import numpy as np
 
-from qiskit.circuit.barrier import Barrier
-from qiskit.circuit.delay import Delay
+from qiskit.circuit import Barrier, Delay, Gate
+from qiskit.circuit.exceptions import CircuitError
 from qiskit.exceptions import QiskitError
 
 
-def _append_circuit(clifford, circuit, qargs=None):
+def _append_circuit(clifford, circuit, qargs=None, recursion_depth=0):
     """Update Clifford inplace by applying a Clifford circuit.
 
     Args:
-        clifford (Clifford): the Clifford to update.
-        circuit (QuantumCircuit): the circuit to apply.
+        clifford (Clifford): The Clifford to update.
+        circuit (QuantumCircuit): The circuit to apply.
         qargs (list or None): The qubits to apply circuit to.
+        recursion_depth (int): The depth of mutual recursion with _append_operation
 
     Returns:
         Clifford: the updated Clifford.
@@ -42,24 +45,26 @@ def _append_circuit(clifford, circuit, qargs=None):
             )
         # Get the integer position of the flat register
         new_qubits = [qargs[circuit.find_bit(bit).index] for bit in instruction.qubits]
-        _append_operation(clifford, instruction.operation, new_qubits)
+        clifford = _append_operation(clifford, instruction.operation, new_qubits, recursion_depth)
     return clifford
 
 
-def _append_operation(clifford, operation, qargs=None):
+def _append_operation(clifford, operation, qargs=None, recursion_depth=0):
     """Update Clifford inplace by applying a Clifford operation.
 
     Args:
-        clifford (Clifford): the Clifford to update.
-        operation (Instruction or str): the operation or composite operation to apply.
+        clifford (Clifford): The Clifford to update.
+        operation (Instruction or Clifford or str): The operation or composite operation to apply.
         qargs (list or None): The qubits to apply operation to.
+        recursion_depth (int): The depth of mutual recursion with _append_circuit
 
     Returns:
         Clifford: the updated Clifford.
 
     Raises:
-        QiskitError: if input operation cannot be decomposed into Clifford operations.
+        QiskitError: if input operation cannot be converted into Clifford operations.
     """
+    # pylint: disable=too-many-return-statements
     if isinstance(operation, (Barrier, Delay)):
         return clifford
 
@@ -91,6 +96,28 @@ def _append_operation(clifford, operation, qargs=None):
             raise QiskitError("Invalid qubits for 2-qubit gate.")
         return _BASIS_2Q[name](clifford, qargs[0], qargs[1])
 
+    # If u gate, check if it is a Clifford, and if so, apply it
+    if isinstance(gate, Gate) and name == "u" and len(qargs) == 1:
+        try:
+            theta, phi, lambd = tuple(_n_half_pis(par) for par in gate.params)
+        except ValueError as err:
+            raise QiskitError("U gate angles must be multiples of pi/2 to be a Clifford") from err
+        if theta == 0:
+            clifford = _append_rz(clifford, qargs[0], lambd + phi)
+        elif theta == 1:
+            clifford = _append_rz(clifford, qargs[0], lambd - 2)
+            clifford = _append_h(clifford, qargs[0])
+            clifford = _append_rz(clifford, qargs[0], phi)
+        elif theta == 2:
+            clifford = _append_rz(clifford, qargs[0], lambd - 1)
+            clifford = _append_x(clifford, qargs[0])
+            clifford = _append_rz(clifford, qargs[0], phi + 1)
+        elif theta == 3:
+            clifford = _append_rz(clifford, qargs[0], lambd)
+            clifford = _append_h(clifford, qargs[0])
+            clifford = _append_rz(clifford, qargs[0], phi + 2)
+        return clifford
+
     # If gate is a Clifford, we can either unroll the gate using the "to_circuit"
     # method, or we can compose the Cliffords directly. Experimentally, for large
     # cliffords the second method is considerably faster.
@@ -103,19 +130,72 @@ def _append_operation(clifford, operation, qargs=None):
         clifford.tableau = composed_clifford.tableau
         return clifford
 
-    # If not a Clifford basis gate we try to unroll the gate and
-    # raise an exception if unrolling reaches a non-Clifford gate.
-    # TODO: We could also check u3 params to see if they
-    # are a single qubit Clifford gate rather than raise an exception.
-    if gate.definition is None:
-        raise QiskitError(f"Cannot apply Instruction: {gate.name}")
-
-    return _append_circuit(clifford, gate.definition, qargs)
+    # If the gate is not directly appendable, we try to unroll the gate with its definition.
+    # This succeeds only if the gate has all-Clifford definition (decomposition).
+    # If fails, we need to restore the clifford that was before attempting to unroll and append.
+    if gate.definition is not None:
+        if recursion_depth > 0:
+            return _append_circuit(clifford, gate.definition, qargs, recursion_depth + 1)
+        else:  # recursion_depth == 0
+            # clifford may be updated in _append_circuit
+            org_clifford = copy.deepcopy(clifford)
+            try:
+                return _append_circuit(clifford, gate.definition, qargs, 1)
+            except (QiskitError, RecursionError):
+                # discard incompletely updated clifford and continue
+                clifford = org_clifford
+
+    # As a final attempt, if the gate is up to 3 qubits,
+    # we try to construct a Clifford to be appended from its matrix representation.
+    if isinstance(gate, Gate) and len(qargs) <= 3:
+        try:
+            matrix = gate.to_matrix()
+            gate_cliff = Clifford.from_matrix(matrix)
+            return _append_operation(clifford, gate_cliff, qargs=qargs)
+        except TypeError as err:
+            raise QiskitError(f"Cannot apply {gate.name} gate with unbounded parameters") from err
+        except CircuitError as err:
+            raise QiskitError(f"Cannot apply {gate.name} gate without to_matrix defined") from err
+        except QiskitError as err:
+            raise QiskitError(f"Cannot apply non-Clifford gate: {gate.name}") from err
+
+    raise QiskitError(f"Cannot apply {gate}")
+
+
+def _n_half_pis(param) -> int:
+    try:
+        param = float(param)
+        epsilon = (abs(param) + 0.5 * 1e-10) % (np.pi / 2)
+        if epsilon > 1e-10:
+            raise ValueError(f"{param} is not to a multiple of pi/2")
+        multiple = int(np.round(param / (np.pi / 2)))
+        return multiple % 4
+    except TypeError as err:
+        raise ValueError(f"{param} is not bounded") from err
 
 
 # ---------------------------------------------------------------------
 # Helper functions for applying basis gates
 # ---------------------------------------------------------------------
+def _append_rz(clifford, qubit, multiple):
+    """Apply an Rz gate to a Clifford.
+
+    Args:
+        clifford (Clifford): a Clifford.
+        qubit (int): gate qubit index.
+        multiple (int): z-rotation angle in a multiple of pi/2
+
+    Returns:
+        Clifford: the updated Clifford.
+    """
+    if multiple % 4 == 1:
+        return _append_s(clifford, qubit)
+    if multiple % 4 == 2:
+        return _append_z(clifford, qubit)
+    if multiple % 4 == 3:
+        return _append_sdg(clifford, qubit)
+
+    return clifford
 
 
 def _append_i(clifford, qubit):

releasenotes/notes/add-clifford-from-matrix-3184822cc559e0b7.yaml

@@ -0,0 +1,9 @@
+---
+features:
+  - |
+    Added :meth:`.Clifford.from_matrix` and :meth:`.Clifford.from_operator` method that
+    creates a ``Clifford`` object from its unitary matrix and operator representation respectively.
+  - |
+    The constructor of :class:`.Clifford` now can take any Clifford gate object up to 3 qubits
+    as long it supports :meth:`to_matrix` method,
+    including parameterized gates such as ``Rz(pi/2)``, which were not convertible before.