diff --git a/cirq-core/cirq/transformers/__init__.py b/cirq-core/cirq/transformers/__init__.py index 5a778095eac..dcb445c3c1b 100644 --- a/cirq-core/cirq/transformers/__init__.py +++ b/cirq-core/cirq/transformers/__init__.py @@ -37,6 +37,7 @@ two_qubit_matrix_to_diagonal_and_cz_operations, two_qubit_matrix_to_ion_operations, two_qubit_matrix_to_sqrt_iswap_operations, + unitary_to_pauli_string, ) from cirq.transformers.heuristic_decompositions import ( diff --git a/cirq-core/cirq/transformers/analytical_decompositions/__init__.py b/cirq-core/cirq/transformers/analytical_decompositions/__init__.py index f25a323a4e0..12f5d3ad0d6 100644 --- a/cirq-core/cirq/transformers/analytical_decompositions/__init__.py +++ b/cirq-core/cirq/transformers/analytical_decompositions/__init__.py @@ -71,3 +71,7 @@ from cirq.transformers.analytical_decompositions.single_to_two_qubit_isometry import ( two_qubit_matrix_to_cz_isometry, ) + +from cirq.transformers.analytical_decompositions.pauli_string_decomposition import ( + unitary_to_pauli_string, +) diff --git a/cirq-core/cirq/transformers/analytical_decompositions/pauli_string_decomposition.py b/cirq-core/cirq/transformers/analytical_decompositions/pauli_string_decomposition.py new file mode 100644 index 00000000000..93b14101288 --- /dev/null +++ b/cirq-core/cirq/transformers/analytical_decompositions/pauli_string_decomposition.py @@ -0,0 +1,130 @@ +# Copyright 2023 The Cirq Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from typing import Optional, Tuple, cast + +import numpy as np +import numpy.typing as npt + +from cirq.ops import DensePauliString +from cirq import protocols + + +def _argmax(V: npt.NDArray) -> Tuple[int, float]: + """Returns a tuple (index of max number, max number).""" + V = (V * V.conj()).real + idx_max = np.argmax(V) + V[idx_max] = 0 + return cast(int, idx_max), np.max(V) + + +def _validate_decomposition(decomposition: DensePauliString, U: npt.NDArray, eps: float) -> bool: + """Returns whether the max absolute value of the elementwise difference is less than eps.""" + got = protocols.unitary(decomposition) + return np.abs(got - U).max() < eps + + +def _fast_walsh_hadamard_transform(V: npt.NDArray) -> None: + """Fast Walsh–Hadamard Transform of an array.""" + m = len(V) + n = m.bit_length() - 1 + for h in [2**i for i in range(n)]: + for i in range(0, m, h * 2): + for j in range(i, i + h): + x = V[j] + y = V[j + h] + V[j] = x + y + V[j + h] = x - y + + +def _conjugate_with_hadamard(U: npt.NDArray) -> npt.NDArray: + """Applies HcUH in O(n4^n) instead of O(8^n).""" + + U = np.copy(U.T) + for i in range(U.shape[1]): + _fast_walsh_hadamard_transform(U[:, i]) + U = U.T + for i in range(U.shape[1]): + _fast_walsh_hadamard_transform(U[:, i]) + return U + + +def unitary_to_pauli_string(U: npt.NDArray, eps: float = 1e-15) -> Optional[DensePauliString]: + """Attempts to find a pauli string (with possible phase) equivalent to U up to eps. + + Based on this answer https://shorturl.at/aA079. + Let x_mask be the index of the maximum number of the first column of U + and z_mask be the index of the maximum number of the first column of H†UH + each of these indicies is n-bits long where U is 2^n x 2^n. + + These two indices/masks encode in binary the indices of the qubits that + have I, X, Y, Z acting on them as follows: + x_mask[i] == 1 and z_mask[i] == 0: X acts on the ith qubit + x_mask[i] == 0 and z_mask[i] == 1: Z acts on the ith qubit + x_mask[i] == 1 and z_mask[i] == 1: Y acts on the ith qubit + x_mask[i] == 0 and z_mask[i] == 0: I acts on the ith qubit + + Args: + U: A square array whose dimension is a power of 2. + eps: numbers smaller than `eps` are considered zero. + + Returns: + A DensePauliString of None. + + Raises: + ValueError: if U is not square with a power of 2 dimension. + """ + + if len(U.shape) != 2 or U.shape[0] != U.shape[1]: + raise ValueError(f'Input has a non-square shape {U}') + n = U.shape[0].bit_length() - 1 + if U.shape[0] != 2**n: + raise ValueError(f'Input dimension {U.shape[0]} isn\'t a power of 2') + + x_msk, second_largest = _argmax(U[:, 0]) + if second_largest > eps: + return None + U_z = _conjugate_with_hadamard(U) + z_msk, second_largest = _argmax(U_z[:, 0]) + if second_largest > eps: + return None + + def select(i): + """Returns the gate that acts on the ith qubit.""" + has_x = (x_msk >> i) & 1 + has_z = (z_msk >> i) & 1 + # The mapping is: + # - has_x and not has_z => X + # - not has_x and has_z => Z + # - has_x and has_z => Y + # - not has_x and not has_z => I + gate_table = ['IX', 'ZY'] + return gate_table[has_z][has_x] + + decomposition = DensePauliString(''.join(select(i) for i in reversed(range(n)))) + + guess = protocols.unitary(decomposition) + if np.abs(guess[x_msk, 0]) < eps: + return None + phase = U[x_msk, 0] / guess[x_msk, 0] + phase /= np.abs(phase) # Make sure |phase| = 1 to avoid rounding issues. + + decomposition = DensePauliString( + ''.join(select(i) for i in reversed(range(n))), coefficient=phase + ) + + if not _validate_decomposition(decomposition, U, eps): + return None + + return decomposition diff --git a/cirq-core/cirq/transformers/analytical_decompositions/pauli_string_decomposition_test.py b/cirq-core/cirq/transformers/analytical_decompositions/pauli_string_decomposition_test.py new file mode 100644 index 00000000000..d1af1d63c2c --- /dev/null +++ b/cirq-core/cirq/transformers/analytical_decompositions/pauli_string_decomposition_test.py @@ -0,0 +1,58 @@ +# Copyright 2023 The Cirq Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from typing import cast +import itertools +import cmath +import pytest + +import numpy as np + +from cirq.ops import DensePauliString, T +from cirq import protocols +from cirq.transformers.analytical_decompositions import unitary_to_pauli_string + + +@pytest.mark.parametrize('phase', [cmath.exp(i * 2 * cmath.pi / 5 * 1j) for i in range(5)]) +@pytest.mark.parametrize( + 'pauli_string', [''.join(p) for p in itertools.product(['', 'I', 'X', 'Y', 'Z'], repeat=4)] +) +def test_unitary_to_pauli_string(pauli_string: str, phase: complex): + want = DensePauliString(pauli_string, coefficient=phase) + got = unitary_to_pauli_string(protocols.unitary(want)) + assert got is not None + assert np.all(want.pauli_mask == got.pauli_mask) + assert np.isclose(cast(np.complex128, want.coefficient), cast(np.complex128, got.coefficient)) + + +def test_unitary_to_pauli_string_non_pauli_input(): + got = unitary_to_pauli_string(protocols.unitary(T)) + assert got is None + + got = unitary_to_pauli_string(np.array([[1, 0], [1, 0]])) + assert got is None + + got = unitary_to_pauli_string(np.array([[1, 1], [0, 2]])) + assert got is None + + got = unitary_to_pauli_string(np.array([[0, 0.5], [1, -1]]), eps=1.1) + assert got is None + + +def test_invalid_input(): + with pytest.raises(ValueError, match='Input has a non-square shape.*'): + _ = unitary_to_pauli_string(np.zeros((2, 3))) + + with pytest.raises(ValueError, match='Input dimension [0-9]* isn\'t a power of 2'): + _ = unitary_to_pauli_string(np.zeros((3, 3)))