From d22f01e67539249ac49fe3053459d1d763241999 Mon Sep 17 00:00:00 2001 From: Eugene Rublenko <16805621+stand-by@users.noreply.github.com> Date: Thu, 11 Jul 2024 21:00:40 -0400 Subject: [PATCH 01/10] One more implementation for sparse pauli strings --- benchmarks/pauli_operations.py | 153 +++++++++++++++++++++++++++++++++ 1 file changed, 153 insertions(+) create mode 100644 benchmarks/pauli_operations.py diff --git a/benchmarks/pauli_operations.py b/benchmarks/pauli_operations.py new file mode 100644 index 0000000..fcb9bd5 --- /dev/null +++ b/benchmarks/pauli_operations.py @@ -0,0 +1,153 @@ +import numpy as np +from dataclasses import dataclass + + +@dataclass +class PauliString: + string: str + weight: float = 1.0 + +@dataclass +class SparseMatrix: + rows: np.ndarray + columns: np.ndarray + values: np.ndarray + +@dataclass +class SparsePauliString: + columns: np.ndarray + values: np.ndarray + weight: float = 1.0 + + +class PauliComposer: + def __init__(self, pauli: PauliString) -> None: + self.pauli = pauli + self.n_qubits = len(pauli.string) + self.n_vals = 1 << self.n_qubits + self.n_ys = pauli.string.count("Y") + + + def __resolve_init_conditions(self) -> None: + first_col = 0 + for p in self.pauli.string: + first_col <<= 1 + if p == "X" or p == "Y": + first_col += 1 + + match self.n_ys % 4: + case 0: + first_val = 1.0 + case 1: + first_val = -1.0j + case 2: + first_val = -1.0 + case 3: + first_val = 1.0j + + return first_col, first_val + + + def sparse_pauli(self) -> SparsePauliString: + cols = np.empty(self.n_vals, dtype=np.int32) + vals = np.empty(self.n_vals, dtype=np.complex128) + cols[0], vals[0] = self.__resolve_init_conditions() + + for l in range(self.n_qubits): + p = self.pauli.string[self.n_qubits - l - 1] + pow_of_two = 1 << l + + new_slice = slice(pow_of_two, 2*pow_of_two) + old_slice = slice(0, pow_of_two) + + match p: + case "I": + cols[new_slice] = cols[old_slice] + pow_of_two + vals[new_slice] = vals[old_slice] + case "X": + cols[new_slice] = cols[old_slice] - pow_of_two + vals[new_slice] = vals[old_slice] + case "Y": + cols[new_slice] = cols[old_slice] - pow_of_two + vals[new_slice] = -vals[old_slice] + case "Z": + cols[new_slice] = cols[old_slice] + pow_of_two + vals[new_slice] = -vals[old_slice] + + return SparsePauliString(weight=self.pauli.weight, columns=cols, values=vals) + + + def sparse_diag_pauli(self) -> SparsePauliString: + assert self.pauli.string.count("X") + self.pauli.string.count("Y") == 0 + + cols = np.arange(self.n_vals, dtype=np.int32) + vals = np.ones(self.n_vals, dtype=np.complex128) + + for l in range(self.n_qubits): + p = self.pauli.string[self.n_qubits - l - 1] + pow_of_two = 1 << l + + new_slice = slice(pow_of_two, 2*pow_of_two) + old_slice = slice(0, pow_of_two) + + match p: + case "I": + vals[new_slice] = vals[old_slice] + case "Z": + vals[new_slice] = -vals[old_slice] + + return SparsePauliString(weight=self.pauli.weight, columns=cols, values=vals) + + + def efficient_sparse_multiply(self, state: np.ndarray) -> np.ndarray: + assert state.ndim == 2 + + cols = np.empty(self.n_vals, dtype=np.int32) + vals = np.empty(self.n_vals, dtype=np.complex128) + cols[0], vals[0] = self.__resolve_init_conditions() + + product = np.empty((self.n_vals, state.shape[1]), dtype=np.complex128) + product[0] = vals[0] * state[cols[0]] + + for l in range(self.n_qubits): + p = self.pauli.string[self.n_qubits - l - 1] + pow_of_two = 1 << l + + new_slice = slice(pow_of_two, 2*pow_of_two) + old_slice = slice(0, pow_of_two) + + match p: + case "I": + cols[new_slice] = cols[old_slice] + pow_of_two + vals[new_slice] = vals[old_slice] + case "X": + cols[new_slice] = cols[old_slice] - pow_of_two + vals[new_slice] = vals[old_slice] + case "Y": + cols[new_slice] = cols[old_slice] - pow_of_two + vals[new_slice] = -vals[old_slice] + case "Z": + cols[new_slice] = cols[old_slice] + pow_of_two + vals[new_slice] = -vals[old_slice] + + vals[new_slice] *= self.pauli.weight + product[new_slice] = vals[new_slice, np.newaxis] * state[cols[new_slice]] + + return product + + +def sparse_pauli_multiply(pauli_string: SparsePauliString, state: np.ndarray) -> np.ndarray: + if state.ndim == 1: + return pauli_string.weight * pauli_string.values * state[pauli_string.columns] + elif state.ndim == 2: + return pauli_string.weight * pauli_string.values[:, np.newaxis] * state[pauli_string.columns] + else: + raise ValueError("state must be a 1D or 2D array") + + +def dense_pauli_string(): + pass + + +def dense_pauli_multiply(): + pass From 4f4520cd0964f1ddd568b9851da16a7b18442714 Mon Sep 17 00:00:00 2001 From: Eugene Rublenko <16805621+stand-by@users.noreply.github.com> Date: Thu, 11 Jul 2024 21:01:38 -0400 Subject: [PATCH 02/10] Bunch of files for future testing and benchmarking --- benchmarks/benchmark_pauli_operations.py | 11 +++++++++++ benchmarks/test_pauli_operations.py | 14 ++++++++++++++ requirements-dev.txt | 4 ++++ 3 files changed, 29 insertions(+) create mode 100644 benchmarks/benchmark_pauli_operations.py create mode 100644 benchmarks/test_pauli_operations.py create mode 100644 requirements-dev.txt diff --git a/benchmarks/benchmark_pauli_operations.py b/benchmarks/benchmark_pauli_operations.py new file mode 100644 index 0000000..9021980 --- /dev/null +++ b/benchmarks/benchmark_pauli_operations.py @@ -0,0 +1,11 @@ +import numpy as np + +from benchmarks.pauli_operations import * + + +def main(): + pass + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/benchmarks/test_pauli_operations.py b/benchmarks/test_pauli_operations.py new file mode 100644 index 0000000..af76cca --- /dev/null +++ b/benchmarks/test_pauli_operations.py @@ -0,0 +1,14 @@ +import pytest +import numpy as np + +from pauli_operations import * + + +def test_sparse_pauli_string(): + assert PauliComposer(PauliString("IZ")).sparse_pauli() == SparsePauliString(1.0, np.array([0,1]), np.array([1.0, 1.0])) + +# TODO test correspondance to strings made from dense mutliplications + + +if __name__ == "__main__": + pytest.main([__file__]) \ No newline at end of file diff --git a/requirements-dev.txt b/requirements-dev.txt new file mode 100644 index 0000000..c813e86 --- /dev/null +++ b/requirements-dev.txt @@ -0,0 +1,4 @@ +numpy +scipy +numba +pytest \ No newline at end of file From a4da0474ef0059c9821765574caab87e5153a179 Mon Sep 17 00:00:00 2001 From: Eugene Rublenko <16805621+stand-by@users.noreply.github.com> Date: Fri, 12 Jul 2024 18:09:49 -0400 Subject: [PATCH 03/10] Missed weight mutliplication in efficient_sparse_multiply --- benchmarks/pauli_operations.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/benchmarks/pauli_operations.py b/benchmarks/pauli_operations.py index fcb9bd5..1ce9df4 100644 --- a/benchmarks/pauli_operations.py +++ b/benchmarks/pauli_operations.py @@ -107,7 +107,7 @@ def efficient_sparse_multiply(self, state: np.ndarray) -> np.ndarray: cols[0], vals[0] = self.__resolve_init_conditions() product = np.empty((self.n_vals, state.shape[1]), dtype=np.complex128) - product[0] = vals[0] * state[cols[0]] + product[0] = self.pauli.weight * vals[0] * state[cols[0]] for l in range(self.n_qubits): p = self.pauli.string[self.n_qubits - l - 1] @@ -130,9 +130,8 @@ def efficient_sparse_multiply(self, state: np.ndarray) -> np.ndarray: cols[new_slice] = cols[old_slice] + pow_of_two vals[new_slice] = -vals[old_slice] - vals[new_slice] *= self.pauli.weight - product[new_slice] = vals[new_slice, np.newaxis] * state[cols[new_slice]] - + product[new_slice] = self.pauli.weight * vals[new_slice, np.newaxis] * state[cols[new_slice]] + return product From db779f0e2bbd2e0e54ab8b9a1b2bafc8b90803cb Mon Sep 17 00:00:00 2001 From: Eugene Rublenko <16805621+stand-by@users.noreply.github.com> Date: Fri, 12 Jul 2024 18:12:36 -0400 Subject: [PATCH 04/10] Methods for dense conversion --- benchmarks/pauli_operations.py | 56 +++++++++++++++++++++------------- 1 file changed, 35 insertions(+), 21 deletions(-) diff --git a/benchmarks/pauli_operations.py b/benchmarks/pauli_operations.py index 1ce9df4..25604fc 100644 --- a/benchmarks/pauli_operations.py +++ b/benchmarks/pauli_operations.py @@ -2,22 +2,53 @@ from dataclasses import dataclass +def pauli_matrices() -> dict: + s0 = np.array([[1,0],[0,1]], dtype=np.complex128) + s1 = np.array([[0,1],[1,0]], dtype=np.complex128) + s2 = np.array([[0,-1j],[1j,0]], dtype=np.complex128) + s3 = np.array([[1,0],[0,-1]], dtype=np.complex128) + return {'I': s0, 'X': s1, 'Y': s2, 'Z': s3, 0: s0, 1: s1, 2: s2, 3: s3} + + @dataclass class PauliString: string: str weight: float = 1.0 + def dense(self) -> np.ndarray: + paulis = pauli_matrices() + matrix = paulis[self.string[-1]] + for p in reversed(self.string[:-1]): + matrix = np.kron(paulis[p], matrix) + return self.weight * matrix + + +# TODO more validation for the shape of inputs @dataclass -class SparseMatrix: - rows: np.ndarray +class SparsePauliString: columns: np.ndarray values: np.ndarray + weight: float = 1.0 + + def multiply(self, state: np.ndarray) -> np.ndarray: + if state.ndim == 1: + return self.weight * self.values * state[self.columns] + elif state.ndim == 2: + return self.weight * self.values[:, np.newaxis] * state[self.columns] + else: + raise ValueError("state must be a 1D or 2D array") + + def dense(self) -> np.ndarray: + matrix = np.zeros((len(self.columns), len(self.columns)), dtype=np.complex128) + matrix[np.arange(len(self.columns)), self.columns] = self.weight * self.values + return matrix + @dataclass -class SparsePauliString: +class SparseMatrix: + rows: np.ndarray columns: np.ndarray values: np.ndarray - weight: float = 1.0 class PauliComposer: @@ -133,20 +164,3 @@ def efficient_sparse_multiply(self, state: np.ndarray) -> np.ndarray: product[new_slice] = self.pauli.weight * vals[new_slice, np.newaxis] * state[cols[new_slice]] return product - - -def sparse_pauli_multiply(pauli_string: SparsePauliString, state: np.ndarray) -> np.ndarray: - if state.ndim == 1: - return pauli_string.weight * pauli_string.values * state[pauli_string.columns] - elif state.ndim == 2: - return pauli_string.weight * pauli_string.values[:, np.newaxis] * state[pauli_string.columns] - else: - raise ValueError("state must be a 1D or 2D array") - - -def dense_pauli_string(): - pass - - -def dense_pauli_multiply(): - pass From 2332f0dfa726d73c83ab8bac197d014801b40696 Mon Sep 17 00:00:00 2001 From: Eugene Rublenko <16805621+stand-by@users.noreply.github.com> Date: Fri, 12 Jul 2024 18:14:43 -0400 Subject: [PATCH 05/10] Tests for python implementation of pauli composer and other aux structures --- benchmarks/test_pauli_operations.py | 176 +++++++++++++++++++++++++++- 1 file changed, 171 insertions(+), 5 deletions(-) diff --git a/benchmarks/test_pauli_operations.py b/benchmarks/test_pauli_operations.py index af76cca..3103823 100644 --- a/benchmarks/test_pauli_operations.py +++ b/benchmarks/test_pauli_operations.py @@ -1,14 +1,180 @@ import pytest import numpy as np +from itertools import permutations, chain -from pauli_operations import * +import pauli_operations +from pauli_operations import PauliString, SparsePauliString, PauliComposer -def test_sparse_pauli_string(): - assert PauliComposer(PauliString("IZ")).sparse_pauli() == SparsePauliString(1.0, np.array([0,1]), np.array([1.0, 1.0])) +@pytest.fixture +def paulis(): + return pauli_operations.pauli_matrices() -# TODO test correspondance to strings made from dense mutliplications + +def test_pauli_strings(paulis): + for p in "IXYZ": + np.testing.assert_array_equal( + PauliString(p).dense(), paulis[p] + ) + + ps = PauliString("III", 0.5) + np.testing.assert_array_equal( + ps.dense(), np.eye(8) * 0.5 + ) + + ps = PauliString(weight=1.0, string="IZ") + np.testing.assert_array_equal( + ps.dense(), + np.kron(paulis["I"], paulis["Z"]) + ) + + ps = PauliString(weight=0.5, string="XYZ") + np.testing.assert_array_equal( + ps.dense(), + np.kron(paulis["X"], np.kron(paulis["Y"], paulis["Z"])) * 0.5 + ) + + ps = SparsePauliString(np.arange(8), np.ones(8), 0.5) + np.testing.assert_array_equal( + ps.dense(), np.eye(8) * 0.5 + ) + m = np.array([[0, 1, 0], + [0, 0, 2], + [3, 0, 0]]) + ps = SparsePauliString(columns=np.array([1, 2, 0]), + values=np.array([1, 2, 3])) + np.testing.assert_array_equal(ps.dense(), m) + + +def test_pauli_composer(paulis): + for p in "IXYZ": + pc = PauliComposer(PauliString(p)) + assert pc.n_qubits == 1 + assert pc.n_vals == 2 + assert pc.n_ys == 0 or p == "Y" + assert pc.n_ys == 1 or p != "Y" + np.testing.assert_array_equal(pc.sparse_pauli().dense(), paulis[p]) + + pc = PauliComposer(PauliString("II", weight=0.2)) + assert pc.n_vals == 4 + np.testing.assert_array_equal(pc.sparse_pauli().dense(), np.eye(4) * 0.2) + np.testing.assert_array_equal(pc.sparse_diag_pauli().dense(), np.eye(4) * 0.2) + + pc = PauliComposer(PauliString("IIII")) + assert pc.n_vals == 16 + np.testing.assert_array_equal(pc.sparse_pauli().dense(), np.eye(16)) + np.testing.assert_array_equal(pc.sparse_diag_pauli().dense(), np.eye(16)) + + pc = PauliComposer(PauliString("II", weight=0.2)) + assert pc.n_vals == 4 + np.testing.assert_array_equal(pc.sparse_pauli().dense(), np.eye(4) * 0.2) + + pc = PauliComposer(PauliString("XXX", weight=1.0)) + assert pc.n_vals == 8 + np.testing.assert_array_equal(pc.sparse_pauli().dense(), np.fliplr(np.eye(8))) + + pc = PauliComposer(PauliString("IY")) + np.testing.assert_array_equal( + pc.sparse_pauli().dense(), + np.block([[paulis['Y'], np.zeros((2,2))], + [np.zeros((2,2)), paulis['Y']]]) + ) + + pc = PauliComposer(PauliString("IZ")) + np.testing.assert_array_equal( + pc.sparse_pauli().dense(), + np.block([[paulis['Z'], np.zeros((2,2))], + [np.zeros((2,2)), paulis['Z']]]) + ) + + +def test_pauli_composer_equivalence(): + rng = np.random.default_rng(321) + + for c in "IXYZ": + w = rng.random() + np.testing.assert_array_equal( + PauliComposer(PauliString(c, w)).sparse_pauli().dense(), + PauliString(c, w).dense() + ) + + for s in permutations("XYZ", 2): + s = ''.join(s) + w = rng.random() + np.testing.assert_array_equal( + PauliComposer(PauliString(s, w)).sparse_pauli().dense(), + PauliString(s, w).dense() + ) + + for s in permutations("IXYZ", 3): + s = ''.join(s) + w = rng.random() + np.testing.assert_array_equal( + PauliComposer(PauliString(s, w)).sparse_pauli().dense(), + PauliString(s, w).dense() + ) + + ixyz = PauliComposer(PauliString("IXYZ")).sparse_pauli().dense() + np.testing.assert_array_equal(ixyz, PauliString("IXYZ").dense()) + + zyxi = PauliComposer(PauliString("ZYXI")).sparse_pauli().dense() + np.testing.assert_array_equal(zyxi, PauliString("ZYXI").dense()) + + assert np.abs(ixyz - zyxi).sum().sum() > 1e-10 + + for s in ["XYIZXYZ", "XXIYYIZZ", "ZIXIZYXX"]: + np.testing.assert_array_equal( + PauliComposer(PauliString(s)).sparse_pauli().dense(), + PauliString(s).dense() + ) + + +def test_sparse_pauli_multiply(): + rng = np.random.default_rng(321) + + for s in chain(list("IXYZ"), list(permutations("IXYZ", 3)), + ["XYIZXYZ", "XXIYYIZZ", "ZIXIZYXX"]): + s = ''.join(s) + n = 2**len(s) + w = rng.random() + psi = rng.random(n) + psi_batch = rng.random((n, 25)) + + np.testing.assert_allclose( + PauliComposer(PauliString(s, w)).sparse_pauli().multiply(psi), + PauliString(s, w).dense().dot(psi), + atol=1e-15 + ) + np.testing.assert_allclose( + PauliComposer(PauliString(s, w)).sparse_pauli().multiply(psi_batch), + PauliString(s, w).dense() @ psi_batch, + atol=1e-15 + ) + + +def test_pauli_composer_multiply(): + rng = np.random.default_rng(321) + + for s in chain(list("IXYZ"), list(permutations("IXYZ", 3)), + ["XYIZXYZ", "XXIYYIZZ", "ZIXIZYXX"]): + s = ''.join(s) + n = 2**len(s) + w = rng.random() + psi = rng.random(n) + psi_batch = rng.random((n, 20)) + + np.testing.assert_allclose( + PauliComposer(PauliString(s, w)).efficient_sparse_multiply( + psi.reshape(-1, 1)).ravel(), + PauliString(s, w).dense().dot(psi), + atol=1e-15 + ) + np.testing.assert_allclose( + PauliComposer(PauliString(s, w)).efficient_sparse_multiply(psi_batch), + PauliString(s, w).dense() @ psi_batch, + atol=1e-15 + ) if __name__ == "__main__": - pytest.main([__file__]) \ No newline at end of file + pytest.main([__file__]) From 07e01af8c628565c1ce35b440e78a2da220e5de5 Mon Sep 17 00:00:00 2001 From: Eugene Rublenko <16805621+stand-by@users.noreply.github.com> Date: Fri, 12 Jul 2024 18:17:08 -0400 Subject: [PATCH 06/10] Remove empty benchmarking file for now --- benchmarks/benchmark_pauli_operations.py | 11 ----------- 1 file changed, 11 deletions(-) delete mode 100644 benchmarks/benchmark_pauli_operations.py diff --git a/benchmarks/benchmark_pauli_operations.py b/benchmarks/benchmark_pauli_operations.py deleted file mode 100644 index 9021980..0000000 --- a/benchmarks/benchmark_pauli_operations.py +++ /dev/null @@ -1,11 +0,0 @@ -import numpy as np - -from benchmarks.pauli_operations import * - - -def main(): - pass - - -if __name__ == "__main__": - main() \ No newline at end of file From 1c1351eb11c214e62385b04b209e2522c38659cc Mon Sep 17 00:00:00 2001 From: Eugene Rublenko <16805621+stand-by@users.noreply.github.com> Date: Mon, 15 Jul 2024 13:12:52 -0400 Subject: [PATCH 07/10] Let ruff do its thing --- benchmarks/pauli_operations.py | 50 ++++++++--------- benchmarks/test_pauli_operations.py | 85 ++++++++++++----------------- 2 files changed, 60 insertions(+), 75 deletions(-) diff --git a/benchmarks/pauli_operations.py b/benchmarks/pauli_operations.py index 25604fc..de5d797 100644 --- a/benchmarks/pauli_operations.py +++ b/benchmarks/pauli_operations.py @@ -3,11 +3,11 @@ def pauli_matrices() -> dict: - s0 = np.array([[1,0],[0,1]], dtype=np.complex128) - s1 = np.array([[0,1],[1,0]], dtype=np.complex128) - s2 = np.array([[0,-1j],[1j,0]], dtype=np.complex128) - s3 = np.array([[1,0],[0,-1]], dtype=np.complex128) - return {'I': s0, 'X': s1, 'Y': s2, 'Z': s3, 0: s0, 1: s1, 2: s2, 3: s3} + s0 = np.array([[1, 0], [0, 1]], dtype=np.complex128) + s1 = np.array([[0, 1], [1, 0]], dtype=np.complex128) + s2 = np.array([[0, -1j], [1j, 0]], dtype=np.complex128) + s3 = np.array([[1, 0], [0, -1]], dtype=np.complex128) + return {"I": s0, "X": s1, "Y": s2, "Z": s3, 0: s0, 1: s1, 2: s2, 3: s3} @dataclass @@ -58,7 +58,6 @@ def __init__(self, pauli: PauliString) -> None: self.n_vals = 1 << self.n_qubits self.n_ys = pauli.string.count("Y") - def __resolve_init_conditions(self) -> None: first_col = 0 for p in self.pauli.string: @@ -75,20 +74,19 @@ def __resolve_init_conditions(self) -> None: first_val = -1.0 case 3: first_val = 1.0j - - return first_col, first_val + return first_col, first_val def sparse_pauli(self) -> SparsePauliString: cols = np.empty(self.n_vals, dtype=np.int32) vals = np.empty(self.n_vals, dtype=np.complex128) cols[0], vals[0] = self.__resolve_init_conditions() - for l in range(self.n_qubits): - p = self.pauli.string[self.n_qubits - l - 1] - pow_of_two = 1 << l - - new_slice = slice(pow_of_two, 2*pow_of_two) + for q in range(self.n_qubits): + p = self.pauli.string[self.n_qubits - q - 1] + pow_of_two = 1 << q + + new_slice = slice(pow_of_two, 2 * pow_of_two) old_slice = slice(0, pow_of_two) match p: @@ -107,18 +105,17 @@ def sparse_pauli(self) -> SparsePauliString: return SparsePauliString(weight=self.pauli.weight, columns=cols, values=vals) - def sparse_diag_pauli(self) -> SparsePauliString: assert self.pauli.string.count("X") + self.pauli.string.count("Y") == 0 cols = np.arange(self.n_vals, dtype=np.int32) vals = np.ones(self.n_vals, dtype=np.complex128) - for l in range(self.n_qubits): - p = self.pauli.string[self.n_qubits - l - 1] - pow_of_two = 1 << l - - new_slice = slice(pow_of_two, 2*pow_of_two) + for q in range(self.n_qubits): + p = self.pauli.string[self.n_qubits - q - 1] + pow_of_two = 1 << q + + new_slice = slice(pow_of_two, 2 * pow_of_two) old_slice = slice(0, pow_of_two) match p: @@ -129,7 +126,6 @@ def sparse_diag_pauli(self) -> SparsePauliString: return SparsePauliString(weight=self.pauli.weight, columns=cols, values=vals) - def efficient_sparse_multiply(self, state: np.ndarray) -> np.ndarray: assert state.ndim == 2 @@ -140,11 +136,11 @@ def efficient_sparse_multiply(self, state: np.ndarray) -> np.ndarray: product = np.empty((self.n_vals, state.shape[1]), dtype=np.complex128) product[0] = self.pauli.weight * vals[0] * state[cols[0]] - for l in range(self.n_qubits): - p = self.pauli.string[self.n_qubits - l - 1] - pow_of_two = 1 << l - - new_slice = slice(pow_of_two, 2*pow_of_two) + for q in range(self.n_qubits): + p = self.pauli.string[self.n_qubits - q - 1] + pow_of_two = 1 << q + + new_slice = slice(pow_of_two, 2 * pow_of_two) old_slice = slice(0, pow_of_two) match p: @@ -161,6 +157,8 @@ def efficient_sparse_multiply(self, state: np.ndarray) -> np.ndarray: cols[new_slice] = cols[old_slice] + pow_of_two vals[new_slice] = -vals[old_slice] - product[new_slice] = self.pauli.weight * vals[new_slice, np.newaxis] * state[cols[new_slice]] + product[new_slice] = ( + self.pauli.weight * vals[new_slice, np.newaxis] * state[cols[new_slice]] + ) return product diff --git a/benchmarks/test_pauli_operations.py b/benchmarks/test_pauli_operations.py index 3103823..d874926 100644 --- a/benchmarks/test_pauli_operations.py +++ b/benchmarks/test_pauli_operations.py @@ -13,36 +13,23 @@ def paulis(): def test_pauli_strings(paulis): for p in "IXYZ": - np.testing.assert_array_equal( - PauliString(p).dense(), paulis[p] - ) + np.testing.assert_array_equal(PauliString(p).dense(), paulis[p]) ps = PauliString("III", 0.5) - np.testing.assert_array_equal( - ps.dense(), np.eye(8) * 0.5 - ) + np.testing.assert_array_equal(ps.dense(), np.eye(8) * 0.5) ps = PauliString(weight=1.0, string="IZ") - np.testing.assert_array_equal( - ps.dense(), - np.kron(paulis["I"], paulis["Z"]) - ) + np.testing.assert_array_equal(ps.dense(), np.kron(paulis["I"], paulis["Z"])) ps = PauliString(weight=0.5, string="XYZ") np.testing.assert_array_equal( - ps.dense(), - np.kron(paulis["X"], np.kron(paulis["Y"], paulis["Z"])) * 0.5 + ps.dense(), np.kron(paulis["X"], np.kron(paulis["Y"], paulis["Z"])) * 0.5 ) ps = SparsePauliString(np.arange(8), np.ones(8), 0.5) - np.testing.assert_array_equal( - ps.dense(), np.eye(8) * 0.5 - ) - m = np.array([[0, 1, 0], - [0, 0, 2], - [3, 0, 0]]) - ps = SparsePauliString(columns=np.array([1, 2, 0]), - values=np.array([1, 2, 3])) + np.testing.assert_array_equal(ps.dense(), np.eye(8) * 0.5) + m = np.array([[0, 1, 0], [0, 0, 2], [3, 0, 0]]) + ps = SparsePauliString(columns=np.array([1, 2, 0]), values=np.array([1, 2, 3])) np.testing.assert_array_equal(ps.dense(), m) @@ -76,15 +63,13 @@ def test_pauli_composer(paulis): pc = PauliComposer(PauliString("IY")) np.testing.assert_array_equal( pc.sparse_pauli().dense(), - np.block([[paulis['Y'], np.zeros((2,2))], - [np.zeros((2,2)), paulis['Y']]]) + np.block([[paulis["Y"], np.zeros((2, 2))], [np.zeros((2, 2)), paulis["Y"]]]), ) - + pc = PauliComposer(PauliString("IZ")) np.testing.assert_array_equal( pc.sparse_pauli().dense(), - np.block([[paulis['Z'], np.zeros((2,2))], - [np.zeros((2,2)), paulis['Z']]]) + np.block([[paulis["Z"], np.zeros((2, 2))], [np.zeros((2, 2)), paulis["Z"]]]), ) @@ -95,23 +80,23 @@ def test_pauli_composer_equivalence(): w = rng.random() np.testing.assert_array_equal( PauliComposer(PauliString(c, w)).sparse_pauli().dense(), - PauliString(c, w).dense() + PauliString(c, w).dense(), ) for s in permutations("XYZ", 2): - s = ''.join(s) + s = "".join(s) w = rng.random() np.testing.assert_array_equal( PauliComposer(PauliString(s, w)).sparse_pauli().dense(), - PauliString(s, w).dense() + PauliString(s, w).dense(), ) for s in permutations("IXYZ", 3): - s = ''.join(s) + s = "".join(s) w = rng.random() np.testing.assert_array_equal( PauliComposer(PauliString(s, w)).sparse_pauli().dense(), - PauliString(s, w).dense() + PauliString(s, w).dense(), ) ixyz = PauliComposer(PauliString("IXYZ")).sparse_pauli().dense() @@ -124,55 +109,57 @@ def test_pauli_composer_equivalence(): for s in ["XYIZXYZ", "XXIYYIZZ", "ZIXIZYXX"]: np.testing.assert_array_equal( - PauliComposer(PauliString(s)).sparse_pauli().dense(), - PauliString(s).dense() + PauliComposer(PauliString(s)).sparse_pauli().dense(), PauliString(s).dense() ) def test_sparse_pauli_multiply(): rng = np.random.default_rng(321) - - for s in chain(list("IXYZ"), list(permutations("IXYZ", 3)), - ["XYIZXYZ", "XXIYYIZZ", "ZIXIZYXX"]): - s = ''.join(s) - n = 2**len(s) + + for s in chain( + list("IXYZ"), list(permutations("IXYZ", 3)), ["XYIZXYZ", "XXIYYIZZ", "ZIXIZYXX"] + ): + s = "".join(s) + n = 2 ** len(s) w = rng.random() psi = rng.random(n) psi_batch = rng.random((n, 25)) - + np.testing.assert_allclose( PauliComposer(PauliString(s, w)).sparse_pauli().multiply(psi), PauliString(s, w).dense().dot(psi), - atol=1e-15 + atol=1e-15, ) np.testing.assert_allclose( PauliComposer(PauliString(s, w)).sparse_pauli().multiply(psi_batch), PauliString(s, w).dense() @ psi_batch, - atol=1e-15 + atol=1e-15, ) def test_pauli_composer_multiply(): rng = np.random.default_rng(321) - - for s in chain(list("IXYZ"), list(permutations("IXYZ", 3)), - ["XYIZXYZ", "XXIYYIZZ", "ZIXIZYXX"]): - s = ''.join(s) - n = 2**len(s) + + for s in chain( + list("IXYZ"), list(permutations("IXYZ", 3)), ["XYIZXYZ", "XXIYYIZZ", "ZIXIZYXX"] + ): + s = "".join(s) + n = 2 ** len(s) w = rng.random() psi = rng.random(n) psi_batch = rng.random((n, 20)) np.testing.assert_allclose( - PauliComposer(PauliString(s, w)).efficient_sparse_multiply( - psi.reshape(-1, 1)).ravel(), + PauliComposer(PauliString(s, w)) + .efficient_sparse_multiply(psi.reshape(-1, 1)) + .ravel(), PauliString(s, w).dense().dot(psi), - atol=1e-15 + atol=1e-15, ) np.testing.assert_allclose( PauliComposer(PauliString(s, w)).efficient_sparse_multiply(psi_batch), PauliString(s, w).dense() @ psi_batch, - atol=1e-15 + atol=1e-15, ) From 060ab5e6e283c9df515b876494431f09974ba1c8 Mon Sep 17 00:00:00 2001 From: Eugene Rublenko <16805621+stand-by@users.noreply.github.com> Date: Mon, 15 Jul 2024 20:42:42 -0400 Subject: [PATCH 08/10] Remove certain things as per comments from James --- benchmarks/pauli_operations.py | 7 ------- requirements-dev.txt | 4 ---- 2 files changed, 11 deletions(-) delete mode 100644 requirements-dev.txt diff --git a/benchmarks/pauli_operations.py b/benchmarks/pauli_operations.py index de5d797..4bccc14 100644 --- a/benchmarks/pauli_operations.py +++ b/benchmarks/pauli_operations.py @@ -44,13 +44,6 @@ def dense(self) -> np.ndarray: return matrix -@dataclass -class SparseMatrix: - rows: np.ndarray - columns: np.ndarray - values: np.ndarray - - class PauliComposer: def __init__(self, pauli: PauliString) -> None: self.pauli = pauli diff --git a/requirements-dev.txt b/requirements-dev.txt deleted file mode 100644 index c813e86..0000000 --- a/requirements-dev.txt +++ /dev/null @@ -1,4 +0,0 @@ -numpy -scipy -numba -pytest \ No newline at end of file From ddc1b30fe86227a71d440f8a44063e373891aade Mon Sep 17 00:00:00 2001 From: Eugene Rublenko <16805621+stand-by@users.noreply.github.com> Date: Wed, 17 Jul 2024 10:02:10 -0400 Subject: [PATCH 09/10] Proper meaning for weight --- benchmarks/pauli_operations.py | 26 +++++----- benchmarks/test_pauli_operations.py | 75 +++++++++++++++-------------- 2 files changed, 50 insertions(+), 51 deletions(-) diff --git a/benchmarks/pauli_operations.py b/benchmarks/pauli_operations.py index 4bccc14..ee1e2a5 100644 --- a/benchmarks/pauli_operations.py +++ b/benchmarks/pauli_operations.py @@ -10,17 +10,18 @@ def pauli_matrices() -> dict: return {"I": s0, "X": s1, "Y": s2, "Z": s3, 0: s0, 1: s1, 2: s2, 3: s3} -@dataclass class PauliString: - string: str - weight: float = 1.0 + def __init__(self, string: str) -> None: + # TODO validate chars in string + self.string = string + self.weight = len(string) - string.count("I") def dense(self) -> np.ndarray: paulis = pauli_matrices() matrix = paulis[self.string[-1]] for p in reversed(self.string[:-1]): matrix = np.kron(paulis[p], matrix) - return self.weight * matrix + return matrix # TODO more validation for the shape of inputs @@ -28,19 +29,18 @@ def dense(self) -> np.ndarray: class SparsePauliString: columns: np.ndarray values: np.ndarray - weight: float = 1.0 def multiply(self, state: np.ndarray) -> np.ndarray: if state.ndim == 1: - return self.weight * self.values * state[self.columns] + return self.values * state[self.columns] elif state.ndim == 2: - return self.weight * self.values[:, np.newaxis] * state[self.columns] + return self.values[:, np.newaxis] * state[self.columns] else: raise ValueError("state must be a 1D or 2D array") def dense(self) -> np.ndarray: matrix = np.zeros((len(self.columns), len(self.columns)), dtype=np.complex128) - matrix[np.arange(len(self.columns)), self.columns] = self.weight * self.values + matrix[np.arange(len(self.columns)), self.columns] = self.values return matrix @@ -96,7 +96,7 @@ def sparse_pauli(self) -> SparsePauliString: cols[new_slice] = cols[old_slice] + pow_of_two vals[new_slice] = -vals[old_slice] - return SparsePauliString(weight=self.pauli.weight, columns=cols, values=vals) + return SparsePauliString(columns=cols, values=vals) def sparse_diag_pauli(self) -> SparsePauliString: assert self.pauli.string.count("X") + self.pauli.string.count("Y") == 0 @@ -117,7 +117,7 @@ def sparse_diag_pauli(self) -> SparsePauliString: case "Z": vals[new_slice] = -vals[old_slice] - return SparsePauliString(weight=self.pauli.weight, columns=cols, values=vals) + return SparsePauliString(columns=cols, values=vals) def efficient_sparse_multiply(self, state: np.ndarray) -> np.ndarray: assert state.ndim == 2 @@ -127,7 +127,7 @@ def efficient_sparse_multiply(self, state: np.ndarray) -> np.ndarray: cols[0], vals[0] = self.__resolve_init_conditions() product = np.empty((self.n_vals, state.shape[1]), dtype=np.complex128) - product[0] = self.pauli.weight * vals[0] * state[cols[0]] + product[0] = vals[0] * state[cols[0]] for q in range(self.n_qubits): p = self.pauli.string[self.n_qubits - q - 1] @@ -150,8 +150,6 @@ def efficient_sparse_multiply(self, state: np.ndarray) -> np.ndarray: cols[new_slice] = cols[old_slice] + pow_of_two vals[new_slice] = -vals[old_slice] - product[new_slice] = ( - self.pauli.weight * vals[new_slice, np.newaxis] * state[cols[new_slice]] - ) + product[new_slice] = vals[new_slice, np.newaxis] * state[cols[new_slice]] return product diff --git a/benchmarks/test_pauli_operations.py b/benchmarks/test_pauli_operations.py index d874926..55ecc81 100644 --- a/benchmarks/test_pauli_operations.py +++ b/benchmarks/test_pauli_operations.py @@ -12,29 +12,37 @@ def paulis(): def test_pauli_strings(paulis): - for p in "IXYZ": + for p in ["I", "X", "Y", "Z"]: + assert PauliString(p).weight == 1 or p == "I" np.testing.assert_array_equal(PauliString(p).dense(), paulis[p]) - ps = PauliString("III", 0.5) - np.testing.assert_array_equal(ps.dense(), np.eye(8) * 0.5) + ps = PauliString("III") + assert ps.weight == 0 + np.testing.assert_array_equal(ps.dense(), np.eye(8)) - ps = PauliString(weight=1.0, string="IZ") + ps = PauliString(string="IZ") + assert ps.weight == 1 np.testing.assert_array_equal(ps.dense(), np.kron(paulis["I"], paulis["Z"])) - ps = PauliString(weight=0.5, string="XYZ") + ps = PauliString(string="XYZ") + assert ps.weight == 3 np.testing.assert_array_equal( - ps.dense(), np.kron(paulis["X"], np.kron(paulis["Y"], paulis["Z"])) * 0.5 + ps.dense(), np.kron(paulis["X"], np.kron(paulis["Y"], paulis["Z"])) ) - ps = SparsePauliString(np.arange(8), np.ones(8), 0.5) - np.testing.assert_array_equal(ps.dense(), np.eye(8) * 0.5) + assert PauliString("XYIZXYZ").weight == 6 + assert PauliString("XXIYYIZZ").weight == 6 + assert PauliString("ZIXIZYXXY").weight == 7 + + ps = SparsePauliString(np.arange(8), np.ones(8)) + np.testing.assert_array_equal(ps.dense(), np.eye(8)) m = np.array([[0, 1, 0], [0, 0, 2], [3, 0, 0]]) ps = SparsePauliString(columns=np.array([1, 2, 0]), values=np.array([1, 2, 3])) np.testing.assert_array_equal(ps.dense(), m) def test_pauli_composer(paulis): - for p in "IXYZ": + for p in ["I", "X", "Y", "Z"]: pc = PauliComposer(PauliString(p)) assert pc.n_qubits == 1 assert pc.n_vals == 2 @@ -42,21 +50,21 @@ def test_pauli_composer(paulis): assert pc.n_ys == 1 or p != "Y" np.testing.assert_array_equal(pc.sparse_pauli().dense(), paulis[p]) - pc = PauliComposer(PauliString("II", weight=0.2)) + pc = PauliComposer(PauliString("II")) assert pc.n_vals == 4 - np.testing.assert_array_equal(pc.sparse_pauli().dense(), np.eye(4) * 0.2) - np.testing.assert_array_equal(pc.sparse_diag_pauli().dense(), np.eye(4) * 0.2) + np.testing.assert_array_equal(pc.sparse_pauli().dense(), np.eye(4)) + np.testing.assert_array_equal(pc.sparse_diag_pauli().dense(), np.eye(4)) pc = PauliComposer(PauliString("IIII")) assert pc.n_vals == 16 np.testing.assert_array_equal(pc.sparse_pauli().dense(), np.eye(16)) np.testing.assert_array_equal(pc.sparse_diag_pauli().dense(), np.eye(16)) - pc = PauliComposer(PauliString("II", weight=0.2)) + pc = PauliComposer(PauliString("II")) assert pc.n_vals == 4 - np.testing.assert_array_equal(pc.sparse_pauli().dense(), np.eye(4) * 0.2) + np.testing.assert_array_equal(pc.sparse_pauli().dense(), np.eye(4)) - pc = PauliComposer(PauliString("XXX", weight=1.0)) + pc = PauliComposer(PauliString("XXX")) assert pc.n_vals == 8 np.testing.assert_array_equal(pc.sparse_pauli().dense(), np.fliplr(np.eye(8))) @@ -74,29 +82,24 @@ def test_pauli_composer(paulis): def test_pauli_composer_equivalence(): - rng = np.random.default_rng(321) - - for c in "IXYZ": - w = rng.random() + for c in ["I", "X", "Y", "Z"]: np.testing.assert_array_equal( - PauliComposer(PauliString(c, w)).sparse_pauli().dense(), - PauliString(c, w).dense(), + PauliComposer(PauliString(c)).sparse_pauli().dense(), + PauliString(c).dense(), ) for s in permutations("XYZ", 2): s = "".join(s) - w = rng.random() np.testing.assert_array_equal( - PauliComposer(PauliString(s, w)).sparse_pauli().dense(), - PauliString(s, w).dense(), + PauliComposer(PauliString(s)).sparse_pauli().dense(), + PauliString(s).dense(), ) for s in permutations("IXYZ", 3): s = "".join(s) - w = rng.random() np.testing.assert_array_equal( - PauliComposer(PauliString(s, w)).sparse_pauli().dense(), - PauliString(s, w).dense(), + PauliComposer(PauliString(s)).sparse_pauli().dense(), + PauliString(s).dense(), ) ixyz = PauliComposer(PauliString("IXYZ")).sparse_pauli().dense() @@ -121,18 +124,17 @@ def test_sparse_pauli_multiply(): ): s = "".join(s) n = 2 ** len(s) - w = rng.random() psi = rng.random(n) psi_batch = rng.random((n, 25)) np.testing.assert_allclose( - PauliComposer(PauliString(s, w)).sparse_pauli().multiply(psi), - PauliString(s, w).dense().dot(psi), + PauliComposer(PauliString(s)).sparse_pauli().multiply(psi), + PauliString(s).dense().dot(psi), atol=1e-15, ) np.testing.assert_allclose( - PauliComposer(PauliString(s, w)).sparse_pauli().multiply(psi_batch), - PauliString(s, w).dense() @ psi_batch, + PauliComposer(PauliString(s)).sparse_pauli().multiply(psi_batch), + PauliString(s).dense() @ psi_batch, atol=1e-15, ) @@ -145,20 +147,19 @@ def test_pauli_composer_multiply(): ): s = "".join(s) n = 2 ** len(s) - w = rng.random() psi = rng.random(n) psi_batch = rng.random((n, 20)) np.testing.assert_allclose( - PauliComposer(PauliString(s, w)) + PauliComposer(PauliString(s)) .efficient_sparse_multiply(psi.reshape(-1, 1)) .ravel(), - PauliString(s, w).dense().dot(psi), + PauliString(s).dense().dot(psi), atol=1e-15, ) np.testing.assert_allclose( - PauliComposer(PauliString(s, w)).efficient_sparse_multiply(psi_batch), - PauliString(s, w).dense() @ psi_batch, + PauliComposer(PauliString(s)).efficient_sparse_multiply(psi_batch), + PauliString(s).dense() @ psi_batch, atol=1e-15, ) From 18e713759d852a78255b04f04b96d1d39f6ebb89 Mon Sep 17 00:00:00 2001 From: Eugene Rublenko <16805621+stand-by@users.noreply.github.com> Date: Wed, 17 Jul 2024 20:39:51 -0400 Subject: [PATCH 10/10] Address comments from James --- benchmarks/pauli_helpers.py | 17 ++ benchmarks/pauli_operations.py | 230 +++++++++++----------------- benchmarks/test_pauli_operations.py | 157 +++++++++---------- 3 files changed, 180 insertions(+), 224 deletions(-) create mode 100644 benchmarks/pauli_helpers.py diff --git a/benchmarks/pauli_helpers.py b/benchmarks/pauli_helpers.py new file mode 100644 index 0000000..27cf0e4 --- /dev/null +++ b/benchmarks/pauli_helpers.py @@ -0,0 +1,17 @@ +import numpy as np + + +def pauli_matrices() -> dict: + s0 = np.array([[1, 0], [0, 1]], dtype=np.complex128) + s1 = np.array([[0, 1], [1, 0]], dtype=np.complex128) + s2 = np.array([[0, -1j], [1j, 0]], dtype=np.complex128) + s3 = np.array([[1, 0], [0, -1]], dtype=np.complex128) + return {"I": s0, "X": s1, "Y": s2, "Z": s3, 0: s0, 1: s1, 2: s2, 3: s3} + + +def naive_pauli_converter(string: str) -> np.ndarray: + paulis = pauli_matrices() + matrix = paulis[string[-1]] + for p in reversed(string[:-1]): + matrix = np.kron(paulis[p], matrix) + return matrix diff --git a/benchmarks/pauli_operations.py b/benchmarks/pauli_operations.py index ee1e2a5..c9db78f 100644 --- a/benchmarks/pauli_operations.py +++ b/benchmarks/pauli_operations.py @@ -1,155 +1,103 @@ import numpy as np -from dataclasses import dataclass - - -def pauli_matrices() -> dict: - s0 = np.array([[1, 0], [0, 1]], dtype=np.complex128) - s1 = np.array([[0, 1], [1, 0]], dtype=np.complex128) - s2 = np.array([[0, -1j], [1j, 0]], dtype=np.complex128) - s3 = np.array([[1, 0], [0, -1]], dtype=np.complex128) - return {"I": s0, "X": s1, "Y": s2, "Z": s3, 0: s0, 1: s1, 2: s2, 3: s3} class PauliString: def __init__(self, string: str) -> None: - # TODO validate chars in string + if not all([c in "IXYZ" for c in string]): + raise ValueError(f"Invalid pauli string {string}") + self.string = string + self.dim = 1 << len(string) self.weight = len(string) - string.count("I") def dense(self) -> np.ndarray: - paulis = pauli_matrices() - matrix = paulis[self.string[-1]] - for p in reversed(self.string[:-1]): - matrix = np.kron(paulis[p], matrix) - return matrix - + columns, values = compose_sparse_pauli(self.string) -# TODO more validation for the shape of inputs -@dataclass -class SparsePauliString: - columns: np.ndarray - values: np.ndarray + matrix = np.zeros((columns.size, values.size), dtype=np.complex128) + matrix[np.arange(columns.size), columns] = values + return matrix def multiply(self, state: np.ndarray) -> np.ndarray: - if state.ndim == 1: - return self.values * state[self.columns] - elif state.ndim == 2: - return self.values[:, np.newaxis] * state[self.columns] - else: - raise ValueError("state must be a 1D or 2D array") + if state.shape[0] != self.dim or state.ndim > 2: + raise ValueError(f"Provided state has inconsistent shape {state.shape}") - def dense(self) -> np.ndarray: - matrix = np.zeros((len(self.columns), len(self.columns)), dtype=np.complex128) - matrix[np.arange(len(self.columns)), self.columns] = self.values - return matrix + columns, values = compose_sparse_pauli(self.string) - -class PauliComposer: - def __init__(self, pauli: PauliString) -> None: - self.pauli = pauli - self.n_qubits = len(pauli.string) - self.n_vals = 1 << self.n_qubits - self.n_ys = pauli.string.count("Y") - - def __resolve_init_conditions(self) -> None: - first_col = 0 - for p in self.pauli.string: - first_col <<= 1 - if p == "X" or p == "Y": - first_col += 1 - - match self.n_ys % 4: - case 0: - first_val = 1.0 - case 1: - first_val = -1.0j - case 2: - first_val = -1.0 - case 3: - first_val = 1.0j - - return first_col, first_val - - def sparse_pauli(self) -> SparsePauliString: - cols = np.empty(self.n_vals, dtype=np.int32) - vals = np.empty(self.n_vals, dtype=np.complex128) - cols[0], vals[0] = self.__resolve_init_conditions() - - for q in range(self.n_qubits): - p = self.pauli.string[self.n_qubits - q - 1] - pow_of_two = 1 << q - - new_slice = slice(pow_of_two, 2 * pow_of_two) - old_slice = slice(0, pow_of_two) - - match p: - case "I": - cols[new_slice] = cols[old_slice] + pow_of_two - vals[new_slice] = vals[old_slice] - case "X": - cols[new_slice] = cols[old_slice] - pow_of_two - vals[new_slice] = vals[old_slice] - case "Y": - cols[new_slice] = cols[old_slice] - pow_of_two - vals[new_slice] = -vals[old_slice] - case "Z": - cols[new_slice] = cols[old_slice] + pow_of_two - vals[new_slice] = -vals[old_slice] - - return SparsePauliString(columns=cols, values=vals) - - def sparse_diag_pauli(self) -> SparsePauliString: - assert self.pauli.string.count("X") + self.pauli.string.count("Y") == 0 - - cols = np.arange(self.n_vals, dtype=np.int32) - vals = np.ones(self.n_vals, dtype=np.complex128) - - for q in range(self.n_qubits): - p = self.pauli.string[self.n_qubits - q - 1] - pow_of_two = 1 << q - - new_slice = slice(pow_of_two, 2 * pow_of_two) - old_slice = slice(0, pow_of_two) - - match p: - case "I": - vals[new_slice] = vals[old_slice] - case "Z": - vals[new_slice] = -vals[old_slice] - - return SparsePauliString(columns=cols, values=vals) - - def efficient_sparse_multiply(self, state: np.ndarray) -> np.ndarray: - assert state.ndim == 2 - - cols = np.empty(self.n_vals, dtype=np.int32) - vals = np.empty(self.n_vals, dtype=np.complex128) - cols[0], vals[0] = self.__resolve_init_conditions() - - product = np.empty((self.n_vals, state.shape[1]), dtype=np.complex128) - product[0] = vals[0] * state[cols[0]] - - for q in range(self.n_qubits): - p = self.pauli.string[self.n_qubits - q - 1] - pow_of_two = 1 << q - - new_slice = slice(pow_of_two, 2 * pow_of_two) - old_slice = slice(0, pow_of_two) - - match p: - case "I": - cols[new_slice] = cols[old_slice] + pow_of_two - vals[new_slice] = vals[old_slice] - case "X": - cols[new_slice] = cols[old_slice] - pow_of_two - vals[new_slice] = vals[old_slice] - case "Y": - cols[new_slice] = cols[old_slice] - pow_of_two - vals[new_slice] = -vals[old_slice] - case "Z": - cols[new_slice] = cols[old_slice] + pow_of_two - vals[new_slice] = -vals[old_slice] - - product[new_slice] = vals[new_slice, np.newaxis] * state[cols[new_slice]] - - return product + if state.ndim == 2: + return values[:, np.newaxis] * state[columns] + else: + return values * state[columns] + + +def compose_sparse_pauli(string: str) -> tuple[np.ndarray, np.ndarray]: + n_qubits = len(string) + n_vals = 1 << n_qubits + n_ys = string.count("Y") + + # initialize cols array with zeros as we need first element to be 0 + cols = np.zeros(n_vals, dtype=np.int32) + vals = np.empty(n_vals, dtype=np.complex128) + + for p in string: + cols[0] <<= 1 + if p == "X" or p == "Y": + cols[0] += 1 + + match n_ys % 4: + case 0: + vals[0] = 1.0 + case 1: + vals[0] = -1.0j + case 2: + vals[0] = -1.0 + case 3: + vals[0] = 1.0j + + for q in range(n_qubits): + p = string[n_qubits - q - 1] + pow_of_two = 1 << q + + new_slice = slice(pow_of_two, 2 * pow_of_two) + old_slice = slice(0, pow_of_two) + + match p: + case "I": + cols[new_slice] = cols[old_slice] + pow_of_two + vals[new_slice] = vals[old_slice] + case "X": + cols[new_slice] = cols[old_slice] - pow_of_two + vals[new_slice] = vals[old_slice] + case "Y": + cols[new_slice] = cols[old_slice] - pow_of_two + vals[new_slice] = -vals[old_slice] + case "Z": + cols[new_slice] = cols[old_slice] + pow_of_two + vals[new_slice] = -vals[old_slice] + + return cols, vals + + +def compose_sparse_diag_pauli(string) -> np.ndarray: + if "X" in string or "Y" in string: + raise ValueError("Pauli string must contain only I and Z characters") + + n_qubits = len(string) + n_vals = 1 << n_qubits + + # initialize vals array with ones as we need first element to be 1 + vals = np.ones(n_vals, dtype=np.complex128) + + for q in range(n_qubits): + p = string[n_qubits - q - 1] + pow_of_two = 1 << q + + new_slice = slice(pow_of_two, 2 * pow_of_two) + old_slice = slice(0, pow_of_two) + + match p: + case "I": + vals[new_slice] = vals[old_slice] + case "Z": + vals[new_slice] = -vals[old_slice] + + return vals diff --git a/benchmarks/test_pauli_operations.py b/benchmarks/test_pauli_operations.py index 55ecc81..df7f300 100644 --- a/benchmarks/test_pauli_operations.py +++ b/benchmarks/test_pauli_operations.py @@ -2,118 +2,134 @@ import numpy as np from itertools import permutations, chain -import pauli_operations -from pauli_operations import PauliString, SparsePauliString, PauliComposer +from pauli_operations import ( + PauliString, + compose_sparse_pauli, + compose_sparse_diag_pauli, +) +from pauli_helpers import pauli_matrices, naive_pauli_converter @pytest.fixture def paulis(): - return pauli_operations.pauli_matrices() + return pauli_matrices() -def test_pauli_strings(paulis): +def test_pauli_string(paulis): for p in ["I", "X", "Y", "Z"]: - assert PauliString(p).weight == 1 or p == "I" - np.testing.assert_array_equal(PauliString(p).dense(), paulis[p]) + ps = PauliString(p) + assert ps.dim == 2 + assert ps.weight == 1 or p == "I" + np.testing.assert_array_equal(ps.dense(), paulis[p]) + np.testing.assert_array_equal(naive_pauli_converter(p), paulis[p]) ps = PauliString("III") + assert ps.dim == 8 assert ps.weight == 0 np.testing.assert_array_equal(ps.dense(), np.eye(8)) + np.testing.assert_array_equal(naive_pauli_converter(ps.string), np.eye(8)) ps = PauliString(string="IZ") + assert ps.dim == 4 assert ps.weight == 1 np.testing.assert_array_equal(ps.dense(), np.kron(paulis["I"], paulis["Z"])) + np.testing.assert_array_equal( + naive_pauli_converter(ps.string), np.kron(paulis["I"], paulis["Z"]) + ) ps = PauliString(string="XYZ") + assert ps.dim == 8 assert ps.weight == 3 np.testing.assert_array_equal( ps.dense(), np.kron(paulis["X"], np.kron(paulis["Y"], paulis["Z"])) ) + np.testing.assert_array_equal( + naive_pauli_converter(ps.string), + np.kron(paulis["X"], np.kron(paulis["Y"], paulis["Z"])), + ) + assert PauliString("XYIZXYZ").dim == 2**7 assert PauliString("XYIZXYZ").weight == 6 + assert PauliString("XXIYYIZZ").dim == 2**8 assert PauliString("XXIYYIZZ").weight == 6 + assert PauliString("ZIXIZYXXY").dim == 2**9 assert PauliString("ZIXIZYXXY").weight == 7 - ps = SparsePauliString(np.arange(8), np.ones(8)) - np.testing.assert_array_equal(ps.dense(), np.eye(8)) - m = np.array([[0, 1, 0], [0, 0, 2], [3, 0, 0]]) - ps = SparsePauliString(columns=np.array([1, 2, 0]), values=np.array([1, 2, 3])) - np.testing.assert_array_equal(ps.dense(), m) +def test_sparse_pauli_composer(paulis): + ps = PauliString("II") + assert ps.dim == 4 + np.testing.assert_array_equal(ps.dense(), np.eye(4)) + np.testing.assert_array_equal(compose_sparse_diag_pauli(ps.string), np.ones(4)) -def test_pauli_composer(paulis): - for p in ["I", "X", "Y", "Z"]: - pc = PauliComposer(PauliString(p)) - assert pc.n_qubits == 1 - assert pc.n_vals == 2 - assert pc.n_ys == 0 or p == "Y" - assert pc.n_ys == 1 or p != "Y" - np.testing.assert_array_equal(pc.sparse_pauli().dense(), paulis[p]) - - pc = PauliComposer(PauliString("II")) - assert pc.n_vals == 4 - np.testing.assert_array_equal(pc.sparse_pauli().dense(), np.eye(4)) - np.testing.assert_array_equal(pc.sparse_diag_pauli().dense(), np.eye(4)) - - pc = PauliComposer(PauliString("IIII")) - assert pc.n_vals == 16 - np.testing.assert_array_equal(pc.sparse_pauli().dense(), np.eye(16)) - np.testing.assert_array_equal(pc.sparse_diag_pauli().dense(), np.eye(16)) - - pc = PauliComposer(PauliString("II")) - assert pc.n_vals == 4 - np.testing.assert_array_equal(pc.sparse_pauli().dense(), np.eye(4)) - - pc = PauliComposer(PauliString("XXX")) - assert pc.n_vals == 8 - np.testing.assert_array_equal(pc.sparse_pauli().dense(), np.fliplr(np.eye(8))) - - pc = PauliComposer(PauliString("IY")) + ps = PauliString("IIII") + assert ps.dim == 16 + np.testing.assert_array_equal(ps.dense(), np.eye(16)) + np.testing.assert_array_equal(compose_sparse_diag_pauli(ps.string), np.ones(16)) + + ps = PauliString("XXX") + assert ps.dim == 8 + np.testing.assert_array_equal(ps.dense(), np.fliplr(np.eye(8))) + + ps = PauliString("IY") np.testing.assert_array_equal( - pc.sparse_pauli().dense(), + ps.dense(), np.block([[paulis["Y"], np.zeros((2, 2))], [np.zeros((2, 2)), paulis["Y"]]]), ) - pc = PauliComposer(PauliString("IZ")) + ps = PauliString("IZ") np.testing.assert_array_equal( - pc.sparse_pauli().dense(), + ps.dense(), np.block([[paulis["Z"], np.zeros((2, 2))], [np.zeros((2, 2)), paulis["Z"]]]), ) + np.testing.assert_array_equal( + compose_sparse_diag_pauli("IZI"), + compose_sparse_pauli("IZI")[1], + ) + + np.testing.assert_array_equal( + compose_sparse_diag_pauli("ZIZI"), + compose_sparse_pauli("ZIZI")[1], + ) -def test_pauli_composer_equivalence(): + np.testing.assert_array_equal( + compose_sparse_diag_pauli("ZZZIII"), + compose_sparse_pauli("ZZZIII")[1], + ) + + +def test_sparse_pauli_composer_equivalence(): for c in ["I", "X", "Y", "Z"]: np.testing.assert_array_equal( - PauliComposer(PauliString(c)).sparse_pauli().dense(), PauliString(c).dense(), + naive_pauli_converter(c), ) for s in permutations("XYZ", 2): s = "".join(s) np.testing.assert_array_equal( - PauliComposer(PauliString(s)).sparse_pauli().dense(), PauliString(s).dense(), + naive_pauli_converter(s), ) for s in permutations("IXYZ", 3): s = "".join(s) np.testing.assert_array_equal( - PauliComposer(PauliString(s)).sparse_pauli().dense(), PauliString(s).dense(), + naive_pauli_converter(s), ) - ixyz = PauliComposer(PauliString("IXYZ")).sparse_pauli().dense() - np.testing.assert_array_equal(ixyz, PauliString("IXYZ").dense()) + ixyz = PauliString("IXYZ").dense() + np.testing.assert_array_equal(ixyz, naive_pauli_converter("IXYZ")) - zyxi = PauliComposer(PauliString("ZYXI")).sparse_pauli().dense() - np.testing.assert_array_equal(zyxi, PauliString("ZYXI").dense()) + zyxi = PauliString("ZYXI").dense() + np.testing.assert_array_equal(zyxi, naive_pauli_converter("ZYXI")) assert np.abs(ixyz - zyxi).sum().sum() > 1e-10 for s in ["XYIZXYZ", "XXIYYIZZ", "ZIXIZYXX"]: - np.testing.assert_array_equal( - PauliComposer(PauliString(s)).sparse_pauli().dense(), PauliString(s).dense() - ) + np.testing.assert_array_equal(PauliString(s).dense(), naive_pauli_converter(s)) def test_sparse_pauli_multiply(): @@ -125,41 +141,16 @@ def test_sparse_pauli_multiply(): s = "".join(s) n = 2 ** len(s) psi = rng.random(n) - psi_batch = rng.random((n, 25)) - - np.testing.assert_allclose( - PauliComposer(PauliString(s)).sparse_pauli().multiply(psi), - PauliString(s).dense().dot(psi), - atol=1e-15, - ) - np.testing.assert_allclose( - PauliComposer(PauliString(s)).sparse_pauli().multiply(psi_batch), - PauliString(s).dense() @ psi_batch, - atol=1e-15, - ) - - -def test_pauli_composer_multiply(): - rng = np.random.default_rng(321) - - for s in chain( - list("IXYZ"), list(permutations("IXYZ", 3)), ["XYIZXYZ", "XXIYYIZZ", "ZIXIZYXX"] - ): - s = "".join(s) - n = 2 ** len(s) - psi = rng.random(n) - psi_batch = rng.random((n, 20)) + psi_batch = rng.random((n, 21)) np.testing.assert_allclose( - PauliComposer(PauliString(s)) - .efficient_sparse_multiply(psi.reshape(-1, 1)) - .ravel(), - PauliString(s).dense().dot(psi), + PauliString(s).multiply(psi), + naive_pauli_converter(s).dot(psi), atol=1e-15, ) np.testing.assert_allclose( - PauliComposer(PauliString(s)).efficient_sparse_multiply(psi_batch), - PauliString(s).dense() @ psi_batch, + PauliString(s).multiply(psi_batch), + naive_pauli_converter(s) @ psi_batch, atol=1e-15, )