From fdbd7516aac9938be5745cd326ca3b5b0d77afe5 Mon Sep 17 00:00:00 2001 From: Jake Lishman Date: Wed, 5 Jan 2022 17:30:14 +0000 Subject: [PATCH 1/3] Speed up Clifford.compose This is a fairly "dumb" optimisation; the algorithm remains exactly the same, but the loops are refactored to do everything meaningful within Numpy vectorised operations. This has a slightly larger memory footprint, than the previous versions in order to vectorise most efficiently, but it is the same time complexity. In practice, however, the speed-ups are very large. On my machine, the timings went (for Cliffords of `n` qubits and made from 10 gates each): n new old 1 0.025 0.057 3 0.146 0.162 10 0.366 1.30 30 1.11 13.5 100 8.51 337 where all timings are in ms. The scaling should not (in theory) be different, but in practice, the faster Numpy vectorisation means larger matrices get bigger speedups. This transforms the 1-qubit case into a lookup table, populated on the first call to `Clifford.compose` with 1-qubit operators. There are 576 possible combinations of 1-qubit Clifford operators, which is small enough to simply look up. Having the special case is an optimisation for randomised benchmarking scripts. (In those cases, doing the lookup manually can save another factor of 2x, because half the time spend in this function is just copying the mutable values from the lookup table.) This also adds a `copy` argument to the constructor of `Clifford` to avoid copying the stabilizer table when it is not necessary. --- .../operators/symplectic/clifford.py | 161 +++++++++++------- ...-compose-performance-96808ba11327e7dd.yaml | 12 ++ 2 files changed, 114 insertions(+), 59 deletions(-) create mode 100644 releasenotes/notes/clifford-compose-performance-96808ba11327e7dd.yaml diff --git a/qiskit/quantum_info/operators/symplectic/clifford.py b/qiskit/quantum_info/operators/symplectic/clifford.py index 167696e8d66c..8acafc929aea 100644 --- a/qiskit/quantum_info/operators/symplectic/clifford.py +++ b/qiskit/quantum_info/operators/symplectic/clifford.py @@ -12,6 +12,7 @@ """ Clifford operator class. """ +import itertools import re import numpy as np @@ -104,18 +105,21 @@ class Clifford(BaseOperator, AdjointMixin, Operation): `arXiv:quant-ph/0406196 `_ """ + _COMPOSE_PHASE_LOOKUP = None + _COMPOSE_1Q_LOOKUP = None + def __array__(self, dtype=None): if dtype: return np.asarray(self.to_matrix(), dtype=dtype) return self.to_matrix() - def __init__(self, data, validate=True): + def __init__(self, data, validate=True, copy=True): """Initialize an operator object.""" - # Initialize from another Clifford by sharing the data + # Initialize from another Clifford if isinstance(data, Clifford): num_qubits = data.num_qubits - self.tableau = data.tableau + self.tableau = data.tableau.copy() if copy else data.tableau # Initialize from ScalarOp as N-qubit identity discarding any global phase elif isinstance(data, ScalarOp): @@ -138,7 +142,7 @@ def __init__(self, data, validate=True): # Initialize StabilizerTable directly from the data else: if isinstance(data, (list, np.ndarray)) and np.asarray(data, dtype=bool).ndim == 2: - data = np.asarray(data, dtype=bool) + data = np.array(data, dtype=bool, copy=copy) if data.shape[0] == data.shape[1]: self.tableau = self._stack_table_phase( data, np.zeros(data.shape[0], dtype=bool) @@ -191,6 +195,9 @@ def __eq__(self, other): """Check if two Clifford tables are equal""" return super().__eq__(other) and (self.tableau == other.tableau).all() + def copy(self): + return type(self)(self, validate=False, copy=True) + # --------------------------------------------------------------------- # Attributes # --------------------------------------------------------------------- @@ -424,7 +431,9 @@ def compose(self, other, qargs=None, front=False): return _append_operation(self.copy(), other, qargs=qargs) if not isinstance(other, Clifford): - other = Clifford(other) + # Not copying is safe since we're going to drop our only reference to `other` at the end + # of the function. + other = Clifford(other, copy=False) # Validate compose dimensions self._op_shape.compose(other._op_shape, qargs, front) @@ -432,64 +441,89 @@ def compose(self, other, qargs=None, front=False): # Pad other with identities if composing on subsystem other = self._pad_with_identity(other, qargs) - if front: - table1 = self - table2 = other - else: - table1 = other - table2 = self - - num_qubits = self.num_qubits - - array1 = table1.symplectic_matrix.astype(int) - phase1 = table1.phase.astype(int) - - array2 = table2.symplectic_matrix.astype(int) - phase2 = table2.phase.astype(int) - - # Update Pauli table - pauli = (array2.dot(array1) % 2).astype(bool) - - # Add phases - phase = np.mod(array2.dot(phase1) + phase2, 2) - - # Correcting for phase due to Pauli multiplication - ifacts = np.zeros(2 * num_qubits, dtype=int) - - for k in range(2 * num_qubits): + left, right = (self, other) if front else (other, self) - row2 = array2[k] - x2 = table2.x[k] - z2 = table2.z[k] - - # Adding a factor of i for each Y in the image of an operator under the - # first operation, since Y=iXZ - - ifacts[k] += np.sum(x2 & z2) - - # Adding factors of i due to qubit-wise Pauli multiplication - - for j in range(num_qubits): - x = 0 - z = 0 - for i in range(2 * num_qubits): - if row2[i]: - x1 = array1[i, j] - z1 = array1[i, j + num_qubits] - if (x | z) & (x1 | z1): - val = np.mod(np.abs(3 * z1 - x1) - np.abs(3 * z - x) - 1, 3) - if val == 0: - ifacts[k] += 1 - elif val == 1: - ifacts[k] -= 1 - x = np.mod(x + x1, 2) - z = np.mod(z + z1, 2) + if self.num_qubits == 1: + return self._compose_1q(left, right) + return self._compose_general(left, right) + @classmethod + def _compose_general(cls, first, second): + # Correcting for phase due to Pauli multiplication. Start with factors of -i from XZ = -iY + # on individual qubits, and then handle multiplication between each qubitwise pair. + ifacts = np.sum(second.x & second.z, axis=1, dtype=int) + + # A lookup table for calculating phases. The indices are + # current_x, current_z, running_x_count, running_z_count + # where all counts taken modulo 2. + lookup = np.zeros((2, 2, 2, 2), dtype=int) + lookup[0, 1, 1, 0] = lookup[1, 0, 1, 1] = lookup[1, 1, 0, 1] = -1 + lookup[0, 1, 1, 1] = lookup[1, 0, 0, 1] = lookup[1, 1, 1, 0] = 1 + + x1, z1 = first.x.astype(np.uint8), first.z.astype(np.uint8) + lookup = cls._compose_lookup() + + # The loop is over 2*n_qubits entries, and the entire loop is cubic in the number of qubits. + for k, row2 in enumerate(second.symplectic_matrix): + x1_select = x1[row2] + z1_select = z1[row2] + x1_accum = np.logical_xor.accumulate(x1_select, axis=0).astype(np.uint8) + z1_accum = np.logical_xor.accumulate(z1_select, axis=0).astype(np.uint8) + indexer = (x1_select[1:], z1_select[1:], x1_accum[:-1], z1_accum[:-1]) + ifacts[k] += np.sum(lookup[indexer]) p = np.mod(ifacts, 4) // 2 - phase = np.mod(phase + p, 2) + phase = ( + (np.matmul(second.symplectic_matrix, first.phase, dtype=int) + second.phase + p) % 2 + ).astype(bool) + data = cls._stack_table_phase( + (np.matmul(second.symplectic_matrix, first.symplectic_matrix, dtype=int) % 2).astype( + bool + ), + phase, + ) + return Clifford(data, validate=False, copy=False) - return Clifford(self._stack_table_phase(pauli, phase), validate=False) + @classmethod + def _compose_1q(cls, first, second): + # 1-qubit composition can be done with a simple lookup table; there are 24 elements in the + # 1q Clifford group, so 576 possible combinations, which is small enough to look up. + if cls._COMPOSE_1Q_LOOKUP is None: + # The valid tables for 1q Cliffords. + tables_1q = np.array( + [ + [[False, True], [True, False]], + [[False, True], [True, True]], + [[True, False], [False, True]], + [[True, False], [True, True]], + [[True, True], [False, True]], + [[True, True], [True, False]], + ] + ) + phases_1q = np.array([[False, False], [False, True], [True, False], [True, True]]) + # Build the lookup table. + cliffords = [ + cls(cls._stack_table_phase(table, phase), validate=False, copy=False) + for table, phase in itertools.product(tables_1q, phases_1q) + ] + cls._COMPOSE_1Q_LOOKUP = { + (cls._hash(left), cls._hash(right)): cls._compose_general(left, right) + for left, right in itertools.product(cliffords, repeat=2) + } + return cls._COMPOSE_1Q_LOOKUP[cls._hash(first), cls._hash(second)].copy() + + @classmethod + def _compose_lookup(cls): + if cls._COMPOSE_PHASE_LOOKUP is None: + # A lookup table for calculating phases. The indices are + # current_x, current_z, running_x_count, running_z_count + # where all counts taken modulo 2. + lookup = np.zeros((2, 2, 2, 2), dtype=int) + lookup[0, 1, 1, 0] = lookup[1, 0, 1, 1] = lookup[1, 1, 0, 1] = -1 + lookup[0, 1, 1, 1] = lookup[1, 0, 0, 1] = lookup[1, 1, 1, 0] = 1 + lookup.setflags(write=False) + cls._COMPOSE_PHASE_LOOKUP = lookup + return cls._COMPOSE_PHASE_LOOKUP # --------------------------------------------------------------------- # Representation conversions @@ -719,6 +753,15 @@ def to_labels(self, array=False, mode="B"): # Internal helper functions # --------------------------------------------------------------------- + def _hash(self): + """Produce a hashable value that is unique for each different Clifford. This should only be + used internally when the classes being hashed are under our control, because classes of this + type are mutable.""" + table = self.table + abits = np.packbits(table.array) + pbits = np.packbits(table.phase) + return (abits.tobytes(), pbits.tobytes()) + @staticmethod def _is_symplectic(mat): """Return True if input is symplectic matrix.""" @@ -768,7 +811,7 @@ def _pad_with_identity(self, clifford, qargs): if qargs is None: return clifford - padded = Clifford(np.eye(2 * self.num_qubits, dtype=bool), validate=False) + padded = Clifford(np.eye(2 * self.num_qubits, dtype=bool), validate=False, copy=False) inds = list(qargs) + [self.num_qubits + i for i in qargs] # Pad Pauli array diff --git a/releasenotes/notes/clifford-compose-performance-96808ba11327e7dd.yaml b/releasenotes/notes/clifford-compose-performance-96808ba11327e7dd.yaml new file mode 100644 index 000000000000..1b7caf6da4dc --- /dev/null +++ b/releasenotes/notes/clifford-compose-performance-96808ba11327e7dd.yaml @@ -0,0 +1,12 @@ +--- +features: + - | + :class:`.Clifford` now takes an optional ``copy`` keyword argument in its + constructor. If set to ``False``, then a :class:`.StabilizerTable` provided + as input will not be copied, but will be used directly. This can have + performance benefits, if the data in the table will never be mutated by any + other means. + - | + The performance of :meth:`.Clifford.compose` has been greatly improved for + all numbers of qubits. For operators of 20 qubits, the speedup is on the + order of 100 times. From 4fb9eaca562099a211636f1409381428bc8b38ee Mon Sep 17 00:00:00 2001 From: Jake Lishman Date: Fri, 2 Dec 2022 15:00:09 +0000 Subject: [PATCH 2/3] Remove duplicated lookup table --- qiskit/quantum_info/operators/symplectic/clifford.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/qiskit/quantum_info/operators/symplectic/clifford.py b/qiskit/quantum_info/operators/symplectic/clifford.py index 8acafc929aea..37c4882a295f 100644 --- a/qiskit/quantum_info/operators/symplectic/clifford.py +++ b/qiskit/quantum_info/operators/symplectic/clifford.py @@ -453,13 +453,6 @@ def _compose_general(cls, first, second): # on individual qubits, and then handle multiplication between each qubitwise pair. ifacts = np.sum(second.x & second.z, axis=1, dtype=int) - # A lookup table for calculating phases. The indices are - # current_x, current_z, running_x_count, running_z_count - # where all counts taken modulo 2. - lookup = np.zeros((2, 2, 2, 2), dtype=int) - lookup[0, 1, 1, 0] = lookup[1, 0, 1, 1] = lookup[1, 1, 0, 1] = -1 - lookup[0, 1, 1, 1] = lookup[1, 0, 0, 1] = lookup[1, 1, 1, 0] = 1 - x1, z1 = first.x.astype(np.uint8), first.z.astype(np.uint8) lookup = cls._compose_lookup() From 6d5024716becb4c37e5196e4668096cd3f58de0d Mon Sep 17 00:00:00 2001 From: Jake Lishman Date: Fri, 2 Dec 2022 15:03:38 +0000 Subject: [PATCH 3/3] Use unified tableau in hash Co-authored-by: Ian Hincks --- qiskit/quantum_info/operators/symplectic/clifford.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/qiskit/quantum_info/operators/symplectic/clifford.py b/qiskit/quantum_info/operators/symplectic/clifford.py index 37c4882a295f..17ec261eb3ec 100644 --- a/qiskit/quantum_info/operators/symplectic/clifford.py +++ b/qiskit/quantum_info/operators/symplectic/clifford.py @@ -750,10 +750,7 @@ def _hash(self): """Produce a hashable value that is unique for each different Clifford. This should only be used internally when the classes being hashed are under our control, because classes of this type are mutable.""" - table = self.table - abits = np.packbits(table.array) - pbits = np.packbits(table.phase) - return (abits.tobytes(), pbits.tobytes()) + return np.packbits(self.tableau).tobytes() @staticmethod def _is_symplectic(mat):