From a4d529ecd1d35d6e0fbf18a49eb4c6ef901a4968 Mon Sep 17 00:00:00 2001 From: Charles MOUSSA Date: Mon, 11 Nov 2024 16:44:26 +0100 Subject: [PATCH 1/8] change shadow algorithm for using only samples --- qadence/measurements/shadow.py | 212 +++++------- .../qadence/test_measurements/test_shadows.py | 323 +++++++----------- 2 files changed, 201 insertions(+), 334 deletions(-) diff --git a/qadence/measurements/shadow.py b/qadence/measurements/shadow.py index a9c25cda..90728051 100644 --- a/qadence/measurements/shadow.py +++ b/qadence/measurements/shadow.py @@ -1,40 +1,33 @@ from __future__ import annotations -from collections import Counter -from functools import reduce - import numpy as np import torch from torch import Tensor from qadence.backend import Backend from qadence.backends.pyqtorch import Backend as PyQBackend -from qadence.blocks import AbstractBlock, chain, kron -from qadence.blocks.block_to_tensor import HMAT, IMAT, SDAGMAT, ZMAT, block_to_tensor +from qadence.blocks import AbstractBlock, kron +from qadence.blocks.block_to_tensor import HMAT, IMAT, SDAGMAT from qadence.blocks.composite import CompositeBlock from qadence.blocks.primitive import PrimitiveBlock from qadence.blocks.utils import get_pauli_blocks, unroll_block_with_scaling from qadence.circuit import QuantumCircuit from qadence.engines.differentiable_backend import DifferentiableBackend +from qadence.measurements.utils import get_qubit_indices_for_op, rotate from qadence.noise import NoiseHandler -from qadence.operations import X, Y, Z +from qadence.operations import I, X, Y, Z from qadence.types import Endianness -from qadence.utils import P0_MATRIX, P1_MATRIX pauli_gates = [X, Y, Z] UNITARY_TENSOR = [ - ZMAT @ HMAT, - SDAGMAT.squeeze(dim=0) @ HMAT, + HMAT, + HMAT @ SDAGMAT, IMAT, ] -def identity(n_qubits: int) -> Tensor: - return torch.eye(2**n_qubits, dtype=torch.complex128) - - def _max_observable_weight(observable: AbstractBlock) -> int: """ Get the maximal weight for the given observable. @@ -88,27 +81,36 @@ def number_of_samples( return N, K -def local_shadow(sample: Counter, unitary_ids: list) -> Tensor: +def nested_operator_indexing( + idx_array: np.ndarray, +) -> list: + """Obtain the list of operators from indices. + + Args: + idx_array (np.ndarray): Indices for obtaining the operators. + + Returns: + list: Map of pauli operators. """ - Compute local shadow by inverting the quantum channel for each projector state. + if idx_array.ndim == 1: + return [pauli_gates[int(ind_pauli)](i) for i, ind_pauli in enumerate(idx_array)] # type: ignore[abstract] + return [nested_operator_indexing(sub_array) for sub_array in idx_array] - See https://arxiv.org/pdf/2002.08953.pdf - Supplementary Material 1 and Eqs. (S17,S44). - Expects a sample bitstring in ILO. +def extract_unitaries(unitary_ids: np.ndarray, n_qubits: int) -> list: + """Sample `shadow_size` pauli strings of `n_qubits`. + + Args: + unitary_ids (np.ndarray): Indices for obtaining the operators. + n_qubits (int): Number of qubits + + Returns: + list: Pauli strings. """ - bitstring = list(sample.keys())[0] - local_density_matrices = [] - for bit, unitary_id in zip(bitstring, unitary_ids): - proj_mat = P0_MATRIX if bit == "0" else P1_MATRIX - unitary_tensor = UNITARY_TENSOR[unitary_id].squeeze(dim=0) - local_density_matrices.append( - 3 * (unitary_tensor.adjoint() @ proj_mat @ unitary_tensor) - identity(1) - ) - if len(local_density_matrices) == 1: - return local_density_matrices[0] - else: - return reduce(torch.kron, local_density_matrices) + unitaries = nested_operator_indexing(unitary_ids) + if n_qubits > 1: + unitaries = [kron(*l) for l in unitaries] + return unitaries def classical_shadow( @@ -119,25 +121,17 @@ def classical_shadow( backend: Backend | DifferentiableBackend = PyQBackend(), noise: NoiseHandler | None = None, endianness: Endianness = Endianness.BIG, -) -> list: - shadow: list = [] - # TODO: Parallelize embarrassingly parallel loop. - for _ in range(shadow_size): - unitary_ids = np.random.randint(0, 3, size=(1, circuit.n_qubits))[0] - random_unitary = [ - pauli_gates[unitary_ids[qubit]](qubit) for qubit in range(circuit.n_qubits) - ] - if len(random_unitary) == 1: - random_unitary_block = random_unitary[0] - else: - random_unitary_block = kron(*random_unitary) - rotated_circuit = QuantumCircuit( - circuit.n_qubits, - chain(circuit.block, random_unitary_block), - ) +) -> tuple[np.ndarray, list[Tensor]]: + unitary_ids = np.random.randint(0, 3, size=(shadow_size, circuit.n_qubits)) + all_unitaries = extract_unitaries(unitary_ids, circuit.n_qubits) + shadow: list = list() + batchsize = state.shape[0] if state else 1 + for i in range(shadow_size): + random_unitary_block = all_unitaries[i] + rotated_circuit = rotate(circuit, (random_unitary_block, 1.0)) # Reverse endianness to get sample bitstrings in ILO. conv_circ = backend.circuit(rotated_circuit) - samples = backend.sample( + batch_samples = backend.sample( circuit=conv_circ, param_values=param_values, n_shots=1, @@ -145,97 +139,62 @@ def classical_shadow( noise=noise, endianness=endianness, ) - batched_shadow = [] - for batch in samples: - batched_shadow.append(local_shadow(sample=batch, unitary_ids=unitary_ids)) - shadow.append(batched_shadow) - - # Reshape the shadow by batches of samples instead of samples of batches. - # FIXME: Improve performance. - return [list(s) for s in zip(*shadow)] - - -def reconstruct_state(shadow: list) -> Tensor: - """Reconstruct the state density matrix for the given shadow.""" - return reduce(torch.add, shadow) / len(shadow) - - -def compute_traces( - qubit_support: tuple, - N: int, - K: int, - shadow: list, - observable: AbstractBlock, - endianness: Endianness = Endianness.BIG, -) -> list: - floor = int(np.floor(N / K)) - traces = [] - # TODO: Parallelize embarrassingly parallel loop. - for k in range(K): - reconstructed_state = reconstruct_state(shadow=shadow[k * floor : (k + 1) * floor]) - # Reshape the observable matrix to fit the density matrix dimensions - # by filling indentites. - # Please note the endianness is also flipped to get results in LE. - # FIXME: Changed below from Little to Big, double-check when Roland is back - # FIXME: Correct these comments. - trace = ( - ( - block_to_tensor( - block=observable, - qubit_support=qubit_support, - endianness=Endianness.BIG, - ).squeeze(dim=0) - @ reconstructed_state - ) - .trace() - .real - ) - traces.append(trace) - return traces + shadow.append(batch_samples) + bitstrings = list() + for b in range(batchsize): + bitstrings.append([list(batch[b].keys())[0] for batch in shadow]) + bitstrings_torch = [ + 1 - 2 * torch.stack([torch.tensor([int(b_i) for b_i in sample]) for sample in batch]) + for batch in bitstrings + ] + return unitary_ids, bitstrings_torch def estimators( - qubit_support: tuple, N: int, K: int, - shadow: list, + unitary_shadow_ids: np.ndarray, + shadow_samples: Tensor, observable: AbstractBlock, - endianness: Endianness = Endianness.BIG, ) -> Tensor: """ - Return estimators (traces of observable times mean density matrix). + Return estimators (aka traces calculated by matching the sampled unitaries to the. - for K equally-sized shadow partitions. + observable and averaging the samples) for K equally-sized shadow partitions. See https://arxiv.org/pdf/2002.08953.pdf Algorithm 1. """ - # If there is no Pauli-Z operator in the observable, - # the sample can't "hit" that measurement. + + obs_qubit_support = observable.qubit_support if isinstance(observable, PrimitiveBlock): - if type(observable) == Z: - traces = compute_traces( - qubit_support=qubit_support, - N=N, - K=K, - shadow=shadow, - observable=observable, - endianness=endianness, - ) - else: - traces = [torch.tensor(0.0)] + if isinstance(observable, I): + return torch.tensor(1.0, dtype=torch.get_default_dtype()) + obs_to_pauli_index = [pauli_gates.index(type(observable))] + elif isinstance(observable, CompositeBlock): - if Z in observable: - traces = compute_traces( - qubit_support=qubit_support, - N=N, - K=K, - shadow=shadow, - observable=observable, - endianness=endianness, - ) + obs_to_pauli_index = [ + pauli_gates.index(type(p)) for p in observable.blocks if not isinstance(p, I) # type: ignore[arg-type] + ] + ind_I = set(get_qubit_indices_for_op((observable, 1.0), I(0))) + obs_qubit_support = tuple([ind for ind in observable.qubit_support if ind not in ind_I]) + + floor = int(np.floor(N / K)) + traces = [] + for k in range(K): + indices_match = np.all( + unitary_shadow_ids[k * floor : (k + 1) * floor, obs_qubit_support] + == obs_to_pauli_index, + axis=1, + ) + if indices_match.sum() > 0: + trace = torch.prod( + shadow_samples[k * floor : (k + 1) * floor][indices_match][:, obs_qubit_support], + axis=-1, + ).sum() / sum(indices_match) + traces.append(trace) else: - traces = [torch.tensor(0.0)] + traces.append(torch.tensor(0.0)) return torch.tensor(traces, dtype=torch.get_default_dtype()) @@ -258,7 +217,7 @@ def estimations( N, K = number_of_samples(observables=observables, accuracy=accuracy, confidence=confidence) if shadow_size is not None: N = shadow_size - shadow = classical_shadow( + unitaries_ids, batch_shadow_samples = classical_shadow( shadow_size=N, circuit=circuit, param_values=param_values, @@ -271,18 +230,17 @@ def estimations( for observable in observables: pauli_decomposition = unroll_block_with_scaling(observable) batch_estimations = [] - for batch in shadow: + for batch in batch_shadow_samples: pauli_term_estimations = [] for pauli_term in pauli_decomposition: # Get the estimators for the current Pauli term. # This is a tensor of size K. estimation = estimators( - qubit_support=circuit.block.qubit_support, N=N, K=K, - shadow=batch, + unitary_shadow_ids=unitaries_ids, + shadow_samples=batch, observable=pauli_term[0], - endianness=endianness, ) # Compute the median of means for the current Pauli term. # Weigh the median by the Pauli term scaling. diff --git a/tests/qadence/test_measurements/test_shadows.py b/tests/qadence/test_measurements/test_shadows.py index d0d8990f..091a6d71 100644 --- a/tests/qadence/test_measurements/test_shadows.py +++ b/tests/qadence/test_measurements/test_shadows.py @@ -1,36 +1,22 @@ from __future__ import annotations -import json -import os -from collections import Counter - import pytest import torch -from torch import Tensor from qadence.backends.api import backend_factory from qadence.blocks.abstract import AbstractBlock -from qadence.blocks.block_to_tensor import IMAT from qadence.blocks.utils import add, chain, kron from qadence.circuit import QuantumCircuit from qadence.constructors import ising_hamiltonian, total_magnetization from qadence.execution import expectation -from qadence.measurements import Measurements from qadence.measurements.shadow import ( - UNITARY_TENSOR, _max_observable_weight, - classical_shadow, estimations, - estimators, - local_shadow, number_of_samples, ) -from qadence.model import QuantumModel from qadence.operations import RX, RY, H, I, X, Y, Z from qadence.parameters import Parameter -from qadence.serialization import deserialize from qadence.types import BackendName, DiffMode -from qadence.utils import P0_MATRIX, P1_MATRIX @pytest.mark.parametrize( @@ -61,86 +47,9 @@ def test_number_of_samples( assert K == exp_samples[1] -@pytest.mark.parametrize( - "sample, unitary_ids, exp_shadow", - [ - ( - Counter({"10": 1}), - [0, 2], - torch.kron( - 3 * (UNITARY_TENSOR[0].adjoint() @ P1_MATRIX @ UNITARY_TENSOR[0]) - IMAT, - 3 * (UNITARY_TENSOR[2].adjoint() @ P0_MATRIX @ UNITARY_TENSOR[2]) - IMAT, - ), - ), - ( - Counter({"0111": 1}), - [2, 0, 2, 2], - torch.kron( - torch.kron( - 3 * (UNITARY_TENSOR[2].adjoint() @ P0_MATRIX @ UNITARY_TENSOR[2]) - IMAT, - 3 * (UNITARY_TENSOR[0].adjoint() @ P1_MATRIX @ UNITARY_TENSOR[0]) - IMAT, - ), - torch.kron( - 3 * (UNITARY_TENSOR[2].adjoint() @ P1_MATRIX @ UNITARY_TENSOR[2]) - IMAT, - 3 * (UNITARY_TENSOR[2].adjoint() @ P1_MATRIX @ UNITARY_TENSOR[2]) - IMAT, - ), - ), - ), - ], -) -def test_local_shadow(sample: Counter, unitary_ids: list, exp_shadow: Tensor) -> None: - shadow = local_shadow(sample=sample, unitary_ids=unitary_ids) - assert torch.allclose(shadow, exp_shadow) - - theta = Parameter("theta") -@pytest.mark.skip(reason="Can't fix the seed for deterministic outputs.") -@pytest.mark.parametrize( - "layer, param_values, exp_shadows", - [ - (X(0) @ X(2), {}, []) - # (kron(RX(0, theta), X(1)), {"theta": torch.tensor([0.5, 1.0, 1.5])}, []) - ], -) -def test_classical_shadow(layer: AbstractBlock, param_values: dict, exp_shadows: list) -> None: - circuit = QuantumCircuit(2, layer) - shadows = classical_shadow( - shadow_size=2, - circuit=circuit, - param_values=param_values, - ) - for shadow, exp_shadow in zip(shadows, exp_shadows): - for batch, exp_batch in zip(shadow, exp_shadow): - assert torch.allclose(batch, exp_batch, atol=1.0e-2) - - -@pytest.mark.parametrize( - "N, K, circuit, param_values, observable, exp_traces", - [ - (2, 1, QuantumCircuit(2, kron(X(0), Z(1))), {}, X(1), torch.tensor([0.0])), - ], -) -def test_estimators( - N: int, - K: int, - circuit: QuantumCircuit, - param_values: dict, - observable: AbstractBlock, - exp_traces: Tensor, -) -> None: - shadows = classical_shadow(shadow_size=N, circuit=circuit, param_values=param_values) - estimated_traces = estimators( - qubit_support=circuit.block.qubit_support, - N=N, - K=K, - shadow=shadows[0], - observable=observable, - ) - assert torch.allclose(estimated_traces, exp_traces) - - @pytest.mark.flaky(max_runs=5) @pytest.mark.parametrize( "circuit, observable, values", @@ -201,119 +110,119 @@ def test_estimations_comparison_exact( } -@pytest.mark.flaky(max_runs=5) -@pytest.mark.parametrize( - "circuit, values, diff_mode", - [ - (QuantumCircuit(2, blocks), values, DiffMode.AD), - (QuantumCircuit(2, blocks), values2, DiffMode.GPSR), - ], -) -def test_estimations_comparison_tomo_forward_pass( - circuit: QuantumCircuit, values: dict, diff_mode: DiffMode -) -> None: - observable = Z(0) ^ circuit.n_qubits - - pyq_backend = backend_factory(BackendName.PYQTORCH, diff_mode=diff_mode) - (conv_circ, conv_obs, embed, params) = pyq_backend.convert(circuit, observable) - pyq_exp_exact = pyq_backend.expectation(conv_circ, conv_obs, embed(params, values)) - model = QuantumModel( - circuit=circuit, - observable=observable, - backend=BackendName.PYQTORCH, - diff_mode=DiffMode.GPSR, - ) - options = {"n_shots": 100000} - estimated_exp_tomo = model.expectation( - values=values, - measurement=Measurements(protocol=Measurements.TOMOGRAPHY, options=options), - ) - new_options = {"accuracy": 0.1, "confidence": 0.1} - estimated_exp_shadow = model.expectation( - values=values, - measurement=Measurements(protocol=Measurements.SHADOW, options=new_options), - ) # N = 54400. - assert torch.allclose(estimated_exp_tomo, pyq_exp_exact, atol=1.0e-2) - assert torch.allclose(estimated_exp_shadow, pyq_exp_exact, atol=0.1) - assert torch.allclose(estimated_exp_shadow, pyq_exp_exact, atol=0.1) - - -@pytest.mark.flaky(max_runs=5) -def test_chemistry_hamiltonian_1() -> None: - from qadence import load - - circuit = load("./tests/test_files/chem_circ.json") - assert isinstance(circuit, QuantumCircuit) - hamiltonian = load("./tests/test_files/chem_ham.json") - assert isinstance(hamiltonian, AbstractBlock) - # Restrict shadow size for faster tests. - kwargs = {"accuracy": 0.1, "confidence": 0.1, "shadow_size": 1000} - param_values = {"theta_0": torch.tensor([1.0])} - - model = QuantumModel( - circuit=circuit, - observable=hamiltonian, - backend=BackendName.PYQTORCH, - diff_mode=DiffMode.GPSR, - ) - exact = model.expectation(values=param_values) - estim = model.expectation( - values=param_values, - measurement=Measurements(protocol=Measurements.SHADOW, options=kwargs), - ) - assert torch.allclose(estim, exact, atol=0.3) - - -@pytest.mark.flaky(max_runs=5) -def test_chemistry_hamiltonian_2() -> None: - from qadence import load - - circuit = load("./tests/test_files/chem_circ.json") - assert isinstance(circuit, QuantumCircuit) - hamiltonian = ising_hamiltonian(2) - assert isinstance(hamiltonian, AbstractBlock) - # Restrict shadow size for faster tests. - kwargs = {"accuracy": 0.1, "confidence": 0.1, "shadow_size": 1000} - param_values = {"theta_0": torch.tensor([1.0])} - - model = QuantumModel( - circuit=circuit, - observable=hamiltonian, - backend=BackendName.PYQTORCH, - diff_mode=DiffMode.GPSR, - ) - exact = model.expectation(values=param_values) - estim = model.expectation( - values=param_values, - measurement=Measurements(protocol=Measurements.SHADOW, options=kwargs), - ) - assert torch.allclose(estim, exact, atol=0.2) - - -def open_chem_obs() -> AbstractBlock: - directory = os.getcwd() - with open(os.path.join(directory, "tests/test_files/h4.json"), "r") as js: - obs = json.loads(js.read()) - return deserialize(obs) # type: ignore[return-value] - - -@pytest.mark.flaky(max_runs=5) -def test_chemistry_hamiltonian_3() -> None: - circuit = QuantumCircuit(4, kron(Z(0), H(1), Z(2), X(3))) - hamiltonian = open_chem_obs() - param_values: dict = dict() - - kwargs = {"accuracy": 0.1, "confidence": 0.1, "shadow_size": 5000} - - model = QuantumModel( - circuit=circuit, - observable=hamiltonian, - backend=BackendName.PYQTORCH, - diff_mode=DiffMode.GPSR, - ) - exact = model.expectation(values=param_values) - estim = model.expectation( - values=param_values, - measurement=Measurements(protocol=Measurements.SHADOW, options=kwargs), - ) - assert torch.allclose(estim, exact, atol=0.3) +# @pytest.mark.flaky(max_runs=5) +# @pytest.mark.parametrize( +# "circuit, values, diff_mode", +# [ +# (QuantumCircuit(2, blocks), values, DiffMode.AD), +# (QuantumCircuit(2, blocks), values2, DiffMode.GPSR), +# ], +# ) +# def test_estimations_comparison_tomo_forward_pass( +# circuit: QuantumCircuit, values: dict, diff_mode: DiffMode +# ) -> None: +# observable = Z(0) ^ circuit.n_qubits + +# pyq_backend = backend_factory(BackendName.PYQTORCH, diff_mode=diff_mode) +# (conv_circ, conv_obs, embed, params) = pyq_backend.convert(circuit, observable) +# pyq_exp_exact = pyq_backend.expectation(conv_circ, conv_obs, embed(params, values)) +# model = QuantumModel( +# circuit=circuit, +# observable=observable, +# backend=BackendName.PYQTORCH, +# diff_mode=DiffMode.GPSR, +# ) +# options = {"n_shots": 100000} +# estimated_exp_tomo = model.expectation( +# values=values, +# measurement=Measurements(protocol=Measurements.TOMOGRAPHY, options=options), +# ) +# new_options = {"accuracy": 0.1, "confidence": 0.1} +# estimated_exp_shadow = model.expectation( +# values=values, +# measurement=Measurements(protocol=Measurements.SHADOW, options=new_options), +# ) # N = 54400. +# assert torch.allclose(estimated_exp_tomo, pyq_exp_exact, atol=1.0e-2) +# assert torch.allclose(estimated_exp_shadow, pyq_exp_exact, atol=0.1) +# assert torch.allclose(estimated_exp_shadow, pyq_exp_exact, atol=0.1) + + +# @pytest.mark.flaky(max_runs=5) +# def test_chemistry_hamiltonian_1() -> None: +# from qadence import load + +# circuit = load("./tests/test_files/chem_circ.json") +# assert isinstance(circuit, QuantumCircuit) +# hamiltonian = load("./tests/test_files/chem_ham.json") +# assert isinstance(hamiltonian, AbstractBlock) +# # Restrict shadow size for faster tests. +# kwargs = {"accuracy": 0.1, "confidence": 0.1, "shadow_size": 1000} +# param_values = {"theta_0": torch.tensor([1.0])} + +# model = QuantumModel( +# circuit=circuit, +# observable=hamiltonian, +# backend=BackendName.PYQTORCH, +# diff_mode=DiffMode.GPSR, +# ) +# exact = model.expectation(values=param_values) +# estim = model.expectation( +# values=param_values, +# measurement=Measurements(protocol=Measurements.SHADOW, options=kwargs), +# ) +# assert torch.allclose(estim, exact, atol=0.3) + + +# @pytest.mark.flaky(max_runs=5) +# def test_chemistry_hamiltonian_2() -> None: +# from qadence import load + +# circuit = load("./tests/test_files/chem_circ.json") +# assert isinstance(circuit, QuantumCircuit) +# hamiltonian = ising_hamiltonian(2) +# assert isinstance(hamiltonian, AbstractBlock) +# # Restrict shadow size for faster tests. +# kwargs = {"accuracy": 0.1, "confidence": 0.1, "shadow_size": 1000} +# param_values = {"theta_0": torch.tensor([1.0])} + +# model = QuantumModel( +# circuit=circuit, +# observable=hamiltonian, +# backend=BackendName.PYQTORCH, +# diff_mode=DiffMode.GPSR, +# ) +# exact = model.expectation(values=param_values) +# estim = model.expectation( +# values=param_values, +# measurement=Measurements(protocol=Measurements.SHADOW, options=kwargs), +# ) +# assert torch.allclose(estim, exact, atol=0.2) + + +# def open_chem_obs() -> AbstractBlock: +# directory = os.getcwd() +# with open(os.path.join(directory, "tests/test_files/h4.json"), "r") as js: +# obs = json.loads(js.read()) +# return deserialize(obs) # type: ignore[return-value] + + +# @pytest.mark.flaky(max_runs=5) +# def test_chemistry_hamiltonian_3() -> None: +# circuit = QuantumCircuit(4, kron(Z(0), H(1), Z(2), X(3))) +# hamiltonian = open_chem_obs() +# param_values: dict = dict() + +# kwargs = {"accuracy": 0.1, "confidence": 0.1, "shadow_size": 5000} + +# model = QuantumModel( +# circuit=circuit, +# observable=hamiltonian, +# backend=BackendName.PYQTORCH, +# diff_mode=DiffMode.GPSR, +# ) +# exact = model.expectation(values=param_values) +# estim = model.expectation( +# values=param_values, +# measurement=Measurements(protocol=Measurements.SHADOW, options=kwargs), +# ) +# assert torch.allclose(estim, exact, atol=0.3) From 86385f88e15bff3be90d33953ec1712fe4955fe3 Mon Sep 17 00:00:00 2001 From: Charles MOUSSA Date: Mon, 11 Nov 2024 16:48:35 +0100 Subject: [PATCH 2/8] fix lint --- qadence/measurements/shadow.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qadence/measurements/shadow.py b/qadence/measurements/shadow.py index 90728051..bf9c9fad 100644 --- a/qadence/measurements/shadow.py +++ b/qadence/measurements/shadow.py @@ -109,7 +109,7 @@ def extract_unitaries(unitary_ids: np.ndarray, n_qubits: int) -> list: """ unitaries = nested_operator_indexing(unitary_ids) if n_qubits > 1: - unitaries = [kron(*l) for l in unitaries] + unitaries = [kron(*list_unitaries) for list_unitaries in unitaries] return unitaries From fb65e0d7ba39ace35b4c4834a66c0b4055543037 Mon Sep 17 00:00:00 2001 From: Charles MOUSSA Date: Mon, 11 Nov 2024 16:56:07 +0100 Subject: [PATCH 3/8] fix if statement on state --- qadence/measurements/shadow.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qadence/measurements/shadow.py b/qadence/measurements/shadow.py index bf9c9fad..740e9d58 100644 --- a/qadence/measurements/shadow.py +++ b/qadence/measurements/shadow.py @@ -125,7 +125,7 @@ def classical_shadow( unitary_ids = np.random.randint(0, 3, size=(shadow_size, circuit.n_qubits)) all_unitaries = extract_unitaries(unitary_ids, circuit.n_qubits) shadow: list = list() - batchsize = state.shape[0] if state else 1 + batchsize = state.shape[0] if state is not None else 1 for i in range(shadow_size): random_unitary_block = all_unitaries[i] rotated_circuit = rotate(circuit, (random_unitary_block, 1.0)) From fc2c432324f702e2388bc3e7a5becd8746a5b743 Mon Sep 17 00:00:00 2001 From: Charles MOUSSA Date: Mon, 11 Nov 2024 17:17:19 +0100 Subject: [PATCH 4/8] calculate batchsize differently --- qadence/measurements/shadow.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/qadence/measurements/shadow.py b/qadence/measurements/shadow.py index 740e9d58..2a7b79e4 100644 --- a/qadence/measurements/shadow.py +++ b/qadence/measurements/shadow.py @@ -125,7 +125,7 @@ def classical_shadow( unitary_ids = np.random.randint(0, 3, size=(shadow_size, circuit.n_qubits)) all_unitaries = extract_unitaries(unitary_ids, circuit.n_qubits) shadow: list = list() - batchsize = state.shape[0] if state is not None else 1 + for i in range(shadow_size): random_unitary_block = all_unitaries[i] rotated_circuit = rotate(circuit, (random_unitary_block, 1.0)) @@ -141,6 +141,7 @@ def classical_shadow( ) shadow.append(batch_samples) bitstrings = list() + batchsize = len(batch_samples) for b in range(batchsize): bitstrings.append([list(batch[b].keys())[0] for batch in shadow]) bitstrings_torch = [ From 9820d946ec99f69af14da3b30fb274115fb69655 Mon Sep 17 00:00:00 2001 From: Charles MOUSSA Date: Mon, 11 Nov 2024 17:53:26 +0100 Subject: [PATCH 5/8] reenable tests chemistry --- .../qadence/test_measurements/test_shadows.py | 238 +++++++++--------- 1 file changed, 122 insertions(+), 116 deletions(-) diff --git a/tests/qadence/test_measurements/test_shadows.py b/tests/qadence/test_measurements/test_shadows.py index 091a6d71..e7021cdd 100644 --- a/tests/qadence/test_measurements/test_shadows.py +++ b/tests/qadence/test_measurements/test_shadows.py @@ -1,5 +1,8 @@ from __future__ import annotations +import json +import os + import pytest import torch @@ -9,13 +12,16 @@ from qadence.circuit import QuantumCircuit from qadence.constructors import ising_hamiltonian, total_magnetization from qadence.execution import expectation +from qadence.measurements import Measurements from qadence.measurements.shadow import ( _max_observable_weight, estimations, number_of_samples, ) +from qadence.model import QuantumModel from qadence.operations import RX, RY, H, I, X, Y, Z from qadence.parameters import Parameter +from qadence.serialization import deserialize from qadence.types import BackendName, DiffMode @@ -110,119 +116,119 @@ def test_estimations_comparison_exact( } -# @pytest.mark.flaky(max_runs=5) -# @pytest.mark.parametrize( -# "circuit, values, diff_mode", -# [ -# (QuantumCircuit(2, blocks), values, DiffMode.AD), -# (QuantumCircuit(2, blocks), values2, DiffMode.GPSR), -# ], -# ) -# def test_estimations_comparison_tomo_forward_pass( -# circuit: QuantumCircuit, values: dict, diff_mode: DiffMode -# ) -> None: -# observable = Z(0) ^ circuit.n_qubits - -# pyq_backend = backend_factory(BackendName.PYQTORCH, diff_mode=diff_mode) -# (conv_circ, conv_obs, embed, params) = pyq_backend.convert(circuit, observable) -# pyq_exp_exact = pyq_backend.expectation(conv_circ, conv_obs, embed(params, values)) -# model = QuantumModel( -# circuit=circuit, -# observable=observable, -# backend=BackendName.PYQTORCH, -# diff_mode=DiffMode.GPSR, -# ) -# options = {"n_shots": 100000} -# estimated_exp_tomo = model.expectation( -# values=values, -# measurement=Measurements(protocol=Measurements.TOMOGRAPHY, options=options), -# ) -# new_options = {"accuracy": 0.1, "confidence": 0.1} -# estimated_exp_shadow = model.expectation( -# values=values, -# measurement=Measurements(protocol=Measurements.SHADOW, options=new_options), -# ) # N = 54400. -# assert torch.allclose(estimated_exp_tomo, pyq_exp_exact, atol=1.0e-2) -# assert torch.allclose(estimated_exp_shadow, pyq_exp_exact, atol=0.1) -# assert torch.allclose(estimated_exp_shadow, pyq_exp_exact, atol=0.1) - - -# @pytest.mark.flaky(max_runs=5) -# def test_chemistry_hamiltonian_1() -> None: -# from qadence import load - -# circuit = load("./tests/test_files/chem_circ.json") -# assert isinstance(circuit, QuantumCircuit) -# hamiltonian = load("./tests/test_files/chem_ham.json") -# assert isinstance(hamiltonian, AbstractBlock) -# # Restrict shadow size for faster tests. -# kwargs = {"accuracy": 0.1, "confidence": 0.1, "shadow_size": 1000} -# param_values = {"theta_0": torch.tensor([1.0])} - -# model = QuantumModel( -# circuit=circuit, -# observable=hamiltonian, -# backend=BackendName.PYQTORCH, -# diff_mode=DiffMode.GPSR, -# ) -# exact = model.expectation(values=param_values) -# estim = model.expectation( -# values=param_values, -# measurement=Measurements(protocol=Measurements.SHADOW, options=kwargs), -# ) -# assert torch.allclose(estim, exact, atol=0.3) - - -# @pytest.mark.flaky(max_runs=5) -# def test_chemistry_hamiltonian_2() -> None: -# from qadence import load - -# circuit = load("./tests/test_files/chem_circ.json") -# assert isinstance(circuit, QuantumCircuit) -# hamiltonian = ising_hamiltonian(2) -# assert isinstance(hamiltonian, AbstractBlock) -# # Restrict shadow size for faster tests. -# kwargs = {"accuracy": 0.1, "confidence": 0.1, "shadow_size": 1000} -# param_values = {"theta_0": torch.tensor([1.0])} - -# model = QuantumModel( -# circuit=circuit, -# observable=hamiltonian, -# backend=BackendName.PYQTORCH, -# diff_mode=DiffMode.GPSR, -# ) -# exact = model.expectation(values=param_values) -# estim = model.expectation( -# values=param_values, -# measurement=Measurements(protocol=Measurements.SHADOW, options=kwargs), -# ) -# assert torch.allclose(estim, exact, atol=0.2) - - -# def open_chem_obs() -> AbstractBlock: -# directory = os.getcwd() -# with open(os.path.join(directory, "tests/test_files/h4.json"), "r") as js: -# obs = json.loads(js.read()) -# return deserialize(obs) # type: ignore[return-value] - - -# @pytest.mark.flaky(max_runs=5) -# def test_chemistry_hamiltonian_3() -> None: -# circuit = QuantumCircuit(4, kron(Z(0), H(1), Z(2), X(3))) -# hamiltonian = open_chem_obs() -# param_values: dict = dict() - -# kwargs = {"accuracy": 0.1, "confidence": 0.1, "shadow_size": 5000} - -# model = QuantumModel( -# circuit=circuit, -# observable=hamiltonian, -# backend=BackendName.PYQTORCH, -# diff_mode=DiffMode.GPSR, -# ) -# exact = model.expectation(values=param_values) -# estim = model.expectation( -# values=param_values, -# measurement=Measurements(protocol=Measurements.SHADOW, options=kwargs), -# ) -# assert torch.allclose(estim, exact, atol=0.3) +@pytest.mark.flaky(max_runs=5) +@pytest.mark.parametrize( + "circuit, values, diff_mode", + [ + (QuantumCircuit(2, blocks), values, DiffMode.AD), + (QuantumCircuit(2, blocks), values2, DiffMode.GPSR), + ], +) +def test_estimations_comparison_tomo_forward_pass( + circuit: QuantumCircuit, values: dict, diff_mode: DiffMode +) -> None: + observable = Z(0) ^ circuit.n_qubits + + pyq_backend = backend_factory(BackendName.PYQTORCH, diff_mode=diff_mode) + (conv_circ, conv_obs, embed, params) = pyq_backend.convert(circuit, observable) + pyq_exp_exact = pyq_backend.expectation(conv_circ, conv_obs, embed(params, values)) + model = QuantumModel( + circuit=circuit, + observable=observable, + backend=BackendName.PYQTORCH, + diff_mode=DiffMode.GPSR, + ) + options = {"n_shots": 100000} + estimated_exp_tomo = model.expectation( + values=values, + measurement=Measurements(protocol=Measurements.TOMOGRAPHY, options=options), + ) + new_options = {"accuracy": 0.1, "confidence": 0.1} + estimated_exp_shadow = model.expectation( + values=values, + measurement=Measurements(protocol=Measurements.SHADOW, options=new_options), + ) # N = 54400. + assert torch.allclose(estimated_exp_tomo, pyq_exp_exact, atol=1.0e-2) + assert torch.allclose(estimated_exp_shadow, pyq_exp_exact, atol=0.1) + assert torch.allclose(estimated_exp_shadow, pyq_exp_exact, atol=0.1) + + +@pytest.mark.flaky(max_runs=5) +def test_chemistry_hamiltonian_1() -> None: + from qadence import load + + circuit = load("./tests/test_files/chem_circ.json") + assert isinstance(circuit, QuantumCircuit) + hamiltonian = load("./tests/test_files/chem_ham.json") + assert isinstance(hamiltonian, AbstractBlock) + # Restrict shadow size for faster tests. + kwargs = {"accuracy": 0.1, "confidence": 0.1, "shadow_size": 1000} + param_values = {"theta_0": torch.tensor([1.0])} + + model = QuantumModel( + circuit=circuit, + observable=hamiltonian, + backend=BackendName.PYQTORCH, + diff_mode=DiffMode.GPSR, + ) + exact = model.expectation(values=param_values) + estim = model.expectation( + values=param_values, + measurement=Measurements(protocol=Measurements.SHADOW, options=kwargs), + ) + assert torch.allclose(estim, exact, atol=0.3) + + +@pytest.mark.flaky(max_runs=5) +def test_chemistry_hamiltonian_2() -> None: + from qadence import load + + circuit = load("./tests/test_files/chem_circ.json") + assert isinstance(circuit, QuantumCircuit) + hamiltonian = ising_hamiltonian(2) + assert isinstance(hamiltonian, AbstractBlock) + # Restrict shadow size for faster tests. + kwargs = {"accuracy": 0.1, "confidence": 0.1, "shadow_size": 1000} + param_values = {"theta_0": torch.tensor([1.0])} + + model = QuantumModel( + circuit=circuit, + observable=hamiltonian, + backend=BackendName.PYQTORCH, + diff_mode=DiffMode.GPSR, + ) + exact = model.expectation(values=param_values) + estim = model.expectation( + values=param_values, + measurement=Measurements(protocol=Measurements.SHADOW, options=kwargs), + ) + assert torch.allclose(estim, exact, atol=0.2) + + +def open_chem_obs() -> AbstractBlock: + directory = os.getcwd() + with open(os.path.join(directory, "tests/test_files/h4.json"), "r") as js: + obs = json.loads(js.read()) + return deserialize(obs) # type: ignore[return-value] + + +@pytest.mark.flaky(max_runs=5) +def test_chemistry_hamiltonian_3() -> None: + circuit = QuantumCircuit(4, kron(Z(0), H(1), Z(2), X(3))) + hamiltonian = open_chem_obs() + param_values: dict = dict() + + kwargs = {"accuracy": 0.1, "confidence": 0.1, "shadow_size": 5000} + + model = QuantumModel( + circuit=circuit, + observable=hamiltonian, + backend=BackendName.PYQTORCH, + diff_mode=DiffMode.GPSR, + ) + exact = model.expectation(values=param_values) + estim = model.expectation( + values=param_values, + measurement=Measurements(protocol=Measurements.SHADOW, options=kwargs), + ) + assert torch.allclose(estim, exact, atol=0.3) From 33671e4be96af58248a0b94022f606a0a3794905 Mon Sep 17 00:00:00 2001 From: Charles MOUSSA Date: Tue, 12 Nov 2024 15:30:46 +0100 Subject: [PATCH 6/8] getting more directly rotations instead of using rotate --- qadence/measurements/shadow.py | 44 +++++++++++++++++++++------------- 1 file changed, 28 insertions(+), 16 deletions(-) diff --git a/qadence/measurements/shadow.py b/qadence/measurements/shadow.py index 2a7b79e4..f1cfaf2f 100644 --- a/qadence/measurements/shadow.py +++ b/qadence/measurements/shadow.py @@ -6,20 +6,24 @@ from qadence.backend import Backend from qadence.backends.pyqtorch import Backend as PyQBackend -from qadence.blocks import AbstractBlock, kron +from qadence.blocks import AbstractBlock, KronBlock, chain, kron from qadence.blocks.block_to_tensor import HMAT, IMAT, SDAGMAT from qadence.blocks.composite import CompositeBlock from qadence.blocks.primitive import PrimitiveBlock from qadence.blocks.utils import get_pauli_blocks, unroll_block_with_scaling from qadence.circuit import QuantumCircuit from qadence.engines.differentiable_backend import DifferentiableBackend -from qadence.measurements.utils import get_qubit_indices_for_op, rotate +from qadence.measurements.utils import get_qubit_indices_for_op from qadence.noise import NoiseHandler -from qadence.operations import I, X, Y, Z +from qadence.operations import H, I, SDagger, X, Y, Z from qadence.types import Endianness pauli_gates = [X, Y, Z] - +pauli_rotations = [ + lambda index: H(index), + lambda index: SDagger(index) * H(index), + lambda index: None, +] UNITARY_TENSOR = [ HMAT, @@ -84,33 +88,37 @@ def number_of_samples( def nested_operator_indexing( idx_array: np.ndarray, ) -> list: - """Obtain the list of operators from indices. + """Obtain the list of rotation operators from indices. Args: idx_array (np.ndarray): Indices for obtaining the operators. Returns: - list: Map of pauli operators. + list: Map of rotations. """ if idx_array.ndim == 1: - return [pauli_gates[int(ind_pauli)](i) for i, ind_pauli in enumerate(idx_array)] # type: ignore[abstract] + return [pauli_rotations[int(ind_pauli)](i) for i, ind_pauli in enumerate(idx_array)] # type: ignore[abstract] return [nested_operator_indexing(sub_array) for sub_array in idx_array] -def extract_unitaries(unitary_ids: np.ndarray, n_qubits: int) -> list: - """Sample `shadow_size` pauli strings of `n_qubits`. +def kron_if_non_empty(list_operations: list) -> KronBlock | None: + filtered_op: list = list(filter(None, list_operations)) + return kron(*filtered_op) if len(filtered_op) > 0 else None + + +def extract_operators(unitary_ids: np.ndarray, n_qubits: int) -> list: + """Sample `shadow_size` rotations of `n_qubits`. Args: unitary_ids (np.ndarray): Indices for obtaining the operators. n_qubits (int): Number of qubits - Returns: list: Pauli strings. """ - unitaries = nested_operator_indexing(unitary_ids) + operations = nested_operator_indexing(unitary_ids) if n_qubits > 1: - unitaries = [kron(*list_unitaries) for list_unitaries in unitaries] - return unitaries + operations = [kron_if_non_empty(ops) for ops in operations] + return operations def classical_shadow( @@ -123,12 +131,16 @@ def classical_shadow( endianness: Endianness = Endianness.BIG, ) -> tuple[np.ndarray, list[Tensor]]: unitary_ids = np.random.randint(0, 3, size=(shadow_size, circuit.n_qubits)) - all_unitaries = extract_unitaries(unitary_ids, circuit.n_qubits) shadow: list = list() + all_rotations = extract_operators(unitary_ids, circuit.n_qubits) for i in range(shadow_size): - random_unitary_block = all_unitaries[i] - rotated_circuit = rotate(circuit, (random_unitary_block, 1.0)) + if all_rotations[i]: + rotated_circuit = QuantumCircuit( + circuit.register, chain(circuit.block, all_rotations[i]) + ) + else: + rotated_circuit = circuit # Reverse endianness to get sample bitstrings in ILO. conv_circ = backend.circuit(rotated_circuit) batch_samples = backend.sample( From 13b1444e1ed6ea448177e80a5f7fd98b45f7ddce Mon Sep 17 00:00:00 2001 From: Charles MOUSSA Date: Tue, 12 Nov 2024 16:16:33 +0100 Subject: [PATCH 7/8] add test X observable --- tests/qadence/test_measurements/test_shadows.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/tests/qadence/test_measurements/test_shadows.py b/tests/qadence/test_measurements/test_shadows.py index e7021cdd..32a26200 100644 --- a/tests/qadence/test_measurements/test_shadows.py +++ b/tests/qadence/test_measurements/test_shadows.py @@ -124,11 +124,13 @@ def test_estimations_comparison_exact( (QuantumCircuit(2, blocks), values2, DiffMode.GPSR), ], ) +@pytest.mark.parametrize("observable", [Z(0) ^ 2, X(1)]) def test_estimations_comparison_tomo_forward_pass( - circuit: QuantumCircuit, values: dict, diff_mode: DiffMode + circuit: QuantumCircuit, + values: dict, + diff_mode: DiffMode, + observable: AbstractBlock, ) -> None: - observable = Z(0) ^ circuit.n_qubits - pyq_backend = backend_factory(BackendName.PYQTORCH, diff_mode=diff_mode) (conv_circ, conv_obs, embed, params) = pyq_backend.convert(circuit, observable) pyq_exp_exact = pyq_backend.expectation(conv_circ, conv_obs, embed(params, values)) From 7fe55d6fe9cd900d73ea88768068c0db67f0de5e Mon Sep 17 00:00:00 2001 From: Charles MOUSSA Date: Tue, 12 Nov 2024 16:31:49 +0100 Subject: [PATCH 8/8] change docstring --- qadence/measurements/shadow.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/qadence/measurements/shadow.py b/qadence/measurements/shadow.py index f1cfaf2f..29978d22 100644 --- a/qadence/measurements/shadow.py +++ b/qadence/measurements/shadow.py @@ -171,9 +171,7 @@ def estimators( observable: AbstractBlock, ) -> Tensor: """ - Return estimators (aka traces calculated by matching the sampled unitaries to the. - - observable and averaging the samples) for K equally-sized shadow partitions. + Return trace estimators from the samples for K equally-sized shadow partitions. See https://arxiv.org/pdf/2002.08953.pdf Algorithm 1.