From 08b56e0065d31d4b05a1274f795f2233e835bc34 Mon Sep 17 00:00:00 2001 From: chMoussa Date: Wed, 20 Nov 2024 08:11:11 +0100 Subject: [PATCH] [Feature] Analog noise (#293) Solves #290 . At the moment, this is a minimal addition in `HamiltonianEvolution` by: - Adding a `noise_operators` argument to `HamiltonianEvolution` - Using `mesolve` or `sesolve` depending on the case. - Incorporating a test - Adding an extra section to the docs. --------- --- docs/noise.md | 40 +++++++++++++++++++++ pyproject.toml | 2 +- pyqtorch/hamiltonians/evolution.py | 56 ++++++++++++++++++++++------- pyqtorch/noise/gates.py | 16 ++++++--- tests/test_analog.py | 58 ++++++++++++++++++++++++++++++ 5 files changed, 154 insertions(+), 18 deletions(-) diff --git a/docs/noise.md b/docs/noise.md index 67e02972..bea18327 100644 --- a/docs/noise.md +++ b/docs/noise.md @@ -123,6 +123,46 @@ def fig_to_html(fig: Figure) -> str: # markdown-exec: hide print(fig_to_html(plt.gcf())) # markdown-exec: hide ``` +## Analog Noise + +Analog noise is made possible by specifying `noise_operators` in `HamiltonianEvolution`: + +```python exec="on" source="material-block" +import torch +from pyqtorch import uniform_state, HamiltonianEvolution +from pyqtorch.matrices import DEFAULT_MATRIX_DTYPE +from pyqtorch.noise import Depolarizing +from pyqtorch.utils import SolverType + +n_qubits = 2 +qubit_targets = list(range(n_qubits)) + +# Random hermitian hamiltonian +matrix = torch.rand(2**n_qubits, 2**n_qubits, dtype=DEFAULT_MATRIX_DTYPE) +hermitian_matrix = matrix + matrix.T.conj() + +time = torch.tensor([1.0]) +time_symbol = "t" +dur_val = torch.rand(1) +list_ops = Depolarizing(0, error_probability=0.1).tensor(2) +list_ops = [op.squeeze() for op in list_ops] +solver = SolverType.DP5_ME +n_steps = 5 + +hamiltonian_evolution = HamiltonianEvolution(hermitian_matrix, time_symbol, qubit_targets, + duration=dur_val, steps=n_steps, + solver=solver, noise_operators=list_ops,) + +# Starting from a uniform state +psi_start = uniform_state(n_qubits) + +# Returns the evolved state +psi_end = hamiltonian_evolution(state = psi_start, values={time_symbol: time}) + +print(psi_end) + + +``` ## Readout errors diff --git a/pyproject.toml b/pyproject.toml index 23daf164..1ac6b19d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,7 +6,7 @@ build-backend = "hatchling.build" name = "pyqtorch" description = "An efficient, large-scale emulator designed for quantum machine learning, seamlessly integrated with a PyTorch backend. Please refer to https://pyqtorch.readthedocs.io/en/latest/ for setup and usage info, along with the full documentation." readme = "README.md" -version = "1.5.2" +version = "1.6.0" requires-python = ">=3.8,<3.13" license = { text = "Apache 2.0" } keywords = ["quantum"] diff --git a/pyqtorch/hamiltonians/evolution.py b/pyqtorch/hamiltonians/evolution.py index 9d46c3af..e41aa040 100644 --- a/pyqtorch/hamiltonians/evolution.py +++ b/pyqtorch/hamiltonians/evolution.py @@ -16,14 +16,17 @@ from pyqtorch.embed import ConcretizedCallable, Embedding from pyqtorch.primitives import Primitive from pyqtorch.quantum_operation import QuantumOperation +from pyqtorch.time_dependent.mesolve import mesolve from pyqtorch.time_dependent.sesolve import sesolve from pyqtorch.utils import ( ATOL, + DensityMatrix, Operator, SolverType, State, StrEnum, _round_operator, + density_mat, expand_operator, finitediff, is_diag, @@ -140,6 +143,11 @@ class HamiltonianEvolution(Sequence): operations: List of operations. cache_length: LRU cache cache_length evolution operators for given set of parameter values. + duration: Total duration for evolving when using a solver. + steps: Number of steps to use when using solver. + solver: Time-dependent Lindblad master equation solver. + noise_operators: List of tensors or Kraus operators adding analog noise + when solving with a Shrodinger equation solver. """ def __init__( @@ -152,6 +160,7 @@ def __init__( steps: int = 100, solver: SolverType = SolverType.DP5_SE, use_sparse: bool = False, + noise_operators: list[Tensor] = list(), ): """Initializes the HamiltonianEvolution. Depending on the generator argument, set the type and set the right generator getter. @@ -162,7 +171,13 @@ def __init__( qubit_support: The qubits the operator acts on. If generator is a quantum operation or sequence of operations, it will be inferred from the generator. - generator_parametric: Whether the generator is parametric or not. + cache_length: LRU cache cache_length evolution operators for given set + of parameter values. + duration: Total duration for evolving when using a solver. + steps: Number of steps to use when using solver. + solver: Time-dependent Lindblad master equation solver. + noise_operators: List of tensors or Kraus operators adding analog noise + when solving with a Shrodinger equation solver. """ self.solver_type = solver @@ -245,6 +260,8 @@ def __init__( self._cache_hamiltonian_evo: dict[str, Tensor] = dict() self.cache_length = cache_length + self.noise_operators: list[Tensor] = noise_operators + @property def generator(self) -> ModuleList: """Returns the operations making the generator. @@ -413,18 +430,33 @@ def Ht(t: torch.Tensor) -> torch.Tensor: .squeeze(2) ) - sol = sesolve( - Ht, - torch.flatten(state, start_dim=0, end_dim=-2), - t_grid, - self.solver_type, - options={"use_sparse": self.use_sparse}, - ) - - # Retrieve the last state of shape (2**n_qubits, batch_size) - state = sol.states[-1] + if len(self.noise_operators) == 0: + sol = sesolve( + Ht, + torch.flatten(state, start_dim=0, end_dim=-2), + t_grid, + self.solver_type, + options={"use_sparse": self.use_sparse}, + ) - return state.reshape([2] * n_qubits + [batch_size]) + # Retrieve the last state of shape (2**n_qubits, batch_size) + # and reshape + state = sol.states[-1].reshape([2] * n_qubits + [batch_size]) + else: + if not isinstance(state, DensityMatrix): + state = density_mat(state) + sol = mesolve( + Ht, + state, + self.noise_operators, + t_grid, + self.solver_type, + options={"use_sparse": self.use_sparse}, + ) + # Retrieve the last density matrix + # and reshape + state = sol.states[-1] + return state def forward( self, diff --git a/pyqtorch/noise/gates.py b/pyqtorch/noise/gates.py index 19a9c7b4..136c0aa9 100644 --- a/pyqtorch/noise/gates.py +++ b/pyqtorch/noise/gates.py @@ -9,7 +9,12 @@ from pyqtorch.apply import apply_operator_dm from pyqtorch.embed import Embedding from pyqtorch.matrices import DEFAULT_MATRIX_DTYPE, IMAT, XMAT, YMAT, ZMAT -from pyqtorch.utils import DensityMatrix, density_mat, qubit_support_as_tuple +from pyqtorch.utils import ( + DensityMatrix, + density_mat, + promote_operator, + qubit_support_as_tuple, +) class Noise(torch.nn.Module): @@ -35,11 +40,12 @@ def extra_repr(self) -> str: def kraus_operators(self) -> list[Tensor]: return [getattr(self, f"kraus_{i}") for i in range(len(self._buffers))] - def tensor( - self, - ) -> list[Tensor]: + def tensor(self, n_qubit_support: int | None = None) -> list[Tensor]: # Since PyQ expects tensor.Size = [2**n_qubits, 2**n_qubits,batch_size]. - return [kraus_op.unsqueeze(2) for kraus_op in self.kraus_operators] + t_ops = [kraus_op.unsqueeze(2) for kraus_op in self.kraus_operators] + if n_qubit_support is None: + return t_ops + return [promote_operator(t, self.target, n_qubit_support) for t in t_ops] def forward( self, diff --git a/tests/test_analog.py b/tests/test_analog.py index 8c0c5527..b21a8c18 100644 --- a/tests/test_analog.py +++ b/tests/test_analog.py @@ -18,10 +18,12 @@ XMAT, ZMAT, ) +from pyqtorch.noise import Depolarizing from pyqtorch.utils import ( ATOL, RTOL, SolverType, + density_mat, is_normalized, operator_kron, overlap, @@ -316,6 +318,7 @@ def test_timedependent( # simulate with time-dependent solver t_points = torch.linspace(0, dur_val[0], n_steps) + psi_solver = pyq.sesolve( torch_hamiltonian, psi_start.reshape(-1, batch_size), @@ -345,6 +348,61 @@ def test_timedependent( assert torch.allclose(psi_solver, psi_hamevo, rtol=RTOL, atol=1.0e-3) +@pytest.mark.parametrize("duration", [torch.rand(1), "duration"]) +@pytest.mark.parametrize( + "batchsize", + [ + 1, + 3, + ], +) +def test_timedependent_with_noise( + tparam: str, + param_y: float, + duration: float, + batchsize: int, + n_steps: int, + torch_hamiltonian: Callable, + hamevo_generator: Sequence, + sin: tuple, + sq: tuple, +) -> None: + + psi_start = density_mat(random_state(2, batchsize)) + dur_val = duration if isinstance(duration, torch.Tensor) else torch.rand(1) + + # simulate with time-dependent solver + t_points = torch.linspace(0, dur_val[0], n_steps) + + # Define jump operators + # Note that we squeeze to remove the batch dimension + list_ops = Depolarizing(0, error_probability=0.1).tensor(2) + list_ops = [op.squeeze() for op in list_ops] + solver = SolverType.DP5_ME + psi_solver = pyq.mesolve( + torch_hamiltonian, psi_start, list_ops, t_points, solver + ).states[-1] + + # simulate with HamiltonianEvolution + embedding = pyq.Embedding( + tparam_name=tparam, + var_to_call={sin[0]: sin[1], sq[0]: sq[1]}, + ) + hamiltonian_evolution = pyq.HamiltonianEvolution( + generator=hamevo_generator, + time=tparam, + duration=duration, + steps=n_steps, + solver=solver, + noise_operators=list_ops, + ) + values = {"y": param_y, "duration": dur_val} + psi_hamevo = hamiltonian_evolution( + state=psi_start, values=values, embedding=embedding + ) + assert torch.allclose(psi_solver, psi_hamevo, rtol=RTOL, atol=1.0e-3) + + @pytest.mark.parametrize("n_qubits", [2, 4, 6]) @pytest.mark.parametrize("batch_size", [1, 2]) def test_hamevo_parametric_gen(n_qubits: int, batch_size: int) -> None: