Skip to content

Commit

Permalink
[Feature] Add AnalogNoise interface (#302)
Browse files Browse the repository at this point in the history
Solves #299 by adding a torch module `AnalogNoise` that keeps the noise
operators, and can expend them to match the qubit support of
`HamiltonianEvolution` when using mesolve. Also differentiate between
digital noise operations.

---------

Co-authored-by: RolandMacDoland <[email protected]>
  • Loading branch information
chMoussa and RolandMacDoland authored Dec 19, 2024
1 parent 7c3e30d commit ef035ec
Show file tree
Hide file tree
Showing 17 changed files with 455 additions and 182 deletions.
28 changes: 20 additions & 8 deletions docs/noise.md
Original file line number Diff line number Diff line change
Expand Up @@ -125,13 +125,13 @@ print(fig_to_html(plt.gcf())) # markdown-exec: hide

## Analog Noise

Analog noise is made possible by specifying `noise_operators` in `HamiltonianEvolution`:
Analog noise is made possible by specifying the `noise` argument in `HamiltonianEvolution` either as a list of tensors defining the jump operators to use when using Schrödinger or Lindblad equation solvers or a `AnalogNoise` instance. An `AnalogNoise` instance can be instantiated by providing a list of tensors, and a `qubit_support` that should be a subset of the qubit support of `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.noise import Depolarizing, AnalogNoise
from pyqtorch.utils import SolverType

n_qubits = 2
Expand All @@ -144,26 +144,38 @@ 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]
noise_ops = Depolarizing(0, error_probability=0.1).tensor(2)
noise_ops = [op.squeeze() for op in noise_ops]
# also can be specified as AnalogNoise
noise_ops = AnalogNoise(noise_ops, qubit_support=(0,1))
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,)
solver=solver, noise=noise_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)
print(psi_end) # markdown-exec: hide
```

There is one predefined `AnalogNoise` available: Depolarizing noise (`AnalogDepolarizing`) defined with jump operators: $L_{0,1,2} = \sqrt{\frac{p}{4}} (X, Y, Z)$. Note we can combine `AnalogNoise` with the `+` operator.



```python exec="on" source="material-block"
from pyqtorch.noise import AnalogDepolarizing
analog_noise = AnalogDepolarizing(error_param=0.1, qubit_support=0) + AnalogDepolarizing(error_param=0.1, qubit_support=1)
# we now have a qubit support acting on qubits 0 and 1
print(analog_noise.qubit_support) # markdown-exec: hide
```


## Readout errors

Another source of noise can be added when performing measurements. This is typically described using confusion matrices of the form:
Expand Down Expand Up @@ -196,6 +208,6 @@ noiseless_expectation = pyq.expectation(circ, state, {"theta": theta}, observabl
readobj = ReadoutNoise(n_qubits, seed=0)
noisycirc = pyq.QuantumCircuit(n_qubits, ops, readobj)
noisy_expectation = pyq.expectation(noisycirc, state, {"theta": theta}, observable=obs, n_shots=1000)
print("Noiseless expectation ", noiseless_expectation.item())
print("Noisy expectation ", noisy_expectation.item())
print(f"Noiseless expectation: {noiseless_expectation.item()}") # markdown-exec: hide
print(f"Noisy expectation: {noisy_expectation.item()}") # markdown-exec: hide
```
2 changes: 1 addition & 1 deletion docs/time_dependent.md
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
## Time-dependent Schrödinger and Lindblad master equation solving

For simulating systems described by time-dependent Hamiltonians `pyqtorch` has a module implementing Schrödinger and Lindblad equation solvers.
For simulating systems described by time-dependent Hamiltonians `pyqtorch` has a module implementing Schrödinger and Lindblad equation solvers.

```python exec="on" source="material-block"
import torch
Expand Down
2 changes: 1 addition & 1 deletion pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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.6.0"
version = "1.7.0"
requires-python = ">=3.8,<3.13"
license = { text = "Apache 2.0" }
keywords = ["quantum"]
Expand Down
30 changes: 15 additions & 15 deletions pyqtorch/hamiltonians/evolution.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,19 +14,17 @@
from pyqtorch.circuit import Sequence
from pyqtorch.composite import Scale
from pyqtorch.embed import ConcretizedCallable, Embedding
from pyqtorch.noise import AnalogNoise
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,
Expand Down Expand Up @@ -160,7 +158,7 @@ def __init__(
steps: int = 100,
solver: SolverType = SolverType.DP5_SE,
use_sparse: bool = False,
noise_operators: list[Tensor] = list(),
noise: list[Tensor] | AnalogNoise | None = None,
):
"""Initializes the HamiltonianEvolution.
Depending on the generator argument, set the type and set the right generator getter.
Expand All @@ -176,7 +174,7 @@ def __init__(
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
noise: List of jump operatoes for noisy simulations or an AnalogNoise
when solving with a Shrodinger equation solver.
"""

Expand Down Expand Up @@ -260,7 +258,14 @@ def __init__(
self._cache_hamiltonian_evo: dict[str, Tensor] = dict()
self.cache_length = cache_length

self.noise_operators: list[Tensor] = noise_operators
if isinstance(noise, list):
noise = AnalogNoise(noise, self.qubit_support)
if noise is not None and set(noise.qubit_support) - set(self.qubit_support):
raise ValueError(
"The noise should be a subset or the same qubit support"
"as HamiltonianEvolution."
)
self.noise = noise

@property
def generator(self) -> ModuleList:
Expand Down Expand Up @@ -430,7 +435,7 @@ def Ht(t: torch.Tensor) -> torch.Tensor:
.squeeze(2)
)

if len(self.noise_operators) == 0:
if self.noise is None:
sol = sesolve(
Ht,
torch.flatten(state, start_dim=0, end_dim=-2),
Expand All @@ -443,19 +448,14 @@ def Ht(t: torch.Tensor) -> torch.Tensor:
# 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(
state,
self.noise_operators,
Ht,
t_grid,
self.solver_type,
options={"use_sparse": self.use_sparse},
full_support=self.qubit_support,
)
# Retrieve the last density matrix
# and reshape
state = sol.states[-1]
return state

def forward(
Expand Down
5 changes: 3 additions & 2 deletions pyqtorch/noise/__init__.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
from __future__ import annotations

from .gates import (
from .analog import AnalogDepolarizing, AnalogNoise
from .digital_gates import (
AmplitudeDamping,
BitFlip,
Depolarizing,
Expand All @@ -10,5 +11,5 @@
PhaseDamping,
PhaseFlip,
)
from .protocol import NoiseProtocol, NoiseType, _repr_noise
from .digital_protocol import DigitalNoiseProtocol, DigitalNoiseType, _repr_noise
from .readout import CorrelatedReadoutNoise, ReadoutNoise, WhiteNoise
200 changes: 200 additions & 0 deletions pyqtorch/noise/analog.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,200 @@
from __future__ import annotations

from math import log, sqrt
from typing import Any, Callable

import torch
from torch import Tensor

from pyqtorch.matrices import XMAT, YMAT, ZMAT
from pyqtorch.qubit_support import Support
from pyqtorch.time_dependent.mesolve import mesolve
from pyqtorch.utils import (
DensityMatrix,
SolverType,
density_mat,
expand_operator,
permute_basis,
)


class AnalogNoise(torch.nn.Module):
"""AnalogNoise is used within `HamiltonianEvolution`
when solving the Schrödinger or a master (Lindblad) equation.
Attributes:
noise_operators (list[Tensor]): The list of jump operators
to use with mesolve. Note you can only define them as
2D tensors, so no batchsize.
qubit_support (int | tuple[int, ...] | Support): The qubits
where the noise_operators are defined over.
"""

def __init__(
self,
noise_operators: list[Tensor],
qubit_support: int | tuple[int, ...] | Support,
) -> None:
super().__init__()
self._qubit_support = (
qubit_support
if isinstance(qubit_support, Support)
else Support(target=qubit_support)
)

for index, tensor in enumerate(noise_operators):
if len(self._qubit_support) < int(log(tensor.shape[1], 2)):
raise ValueError(
"Tensor dimensions should match the length of the qubit support."
f"Tensor {tensor} has incompatible dimensions {int(log(tensor.shape[1], 2))}."
)
self.register_buffer(f"noise_operators_{index}", tensor)

self._device: torch.device = noise_operators[0].device
self._dtype: torch.dtype = noise_operators[0].dtype

@property
def noise_operators(self) -> list[Tensor]:
"""Get noise_operators.
Returns:
list[Tensor]: Noise operators.
"""
return [
getattr(self, f"noise_operators_{i}") for i in range(len(self._buffers))
]

@property
def qubit_support(self) -> tuple[int, ...]:
"""Getter qubit_support.
Returns:
Support: Tuple of sorted qubits.
"""
return self._qubit_support.sorted_qubits

def extra_repr(self) -> str:
return f"noise_operators: {self.noise_operators}"

def _noise_operators(
self,
full_support: tuple[int, ...] | None = None,
) -> list[Tensor]:
"""Obtain noise operators expended with full support.
Args:
full_support (tuple[int, ...] | None, optional): _description_. Defaults to None.
Returns:
list[Tensor]: Noise operators defined over `full_support`.
"""
list_ops = self.noise_operators
if self._qubit_support.qubits != self.qubit_support:
list_ops = [
permute_basis(
blockmat.unsqueeze(-1), self._qubit_support.qubits, inv=True
).squeeze(-1)
for blockmat in list_ops
]
if full_support is None:
return list_ops
else:
return [
expand_operator(
blockmat.unsqueeze(-1), self.qubit_support, full_support
).squeeze(-1)
for blockmat in list_ops
]

def to(self, *args: Any, **kwargs: Any) -> AnalogNoise:
"""Do device or dtype conversions.
Returns:
AnalogNoise: Converted instance.
"""
super().to(*args, **kwargs)
self._device = self.noise_operators_0.device
self._dtype = self.noise_operators_0.dtype
return self

def forward(
self,
state: Tensor,
H: Callable[..., Any],
tsave: list | Tensor,
solver: SolverType,
options: dict[str, Any] | None = None,
full_support: tuple[int, ...] | None = None,
) -> Tensor:
"""Obtain the output density matrix by solving a Shrodinger equation.
This uses the `mesolve` function.
Args:
state (Tensor): Input state or density matrix.
H (Callable[..., Any]): Time-dependent hamiltonian fonction.
tsave (list | Tensor): Tensor containing simulation time instants.
solver (SolverType): Name of the solver to use.
options (dict[str, Any] | None, optional): Additional options passed to the solver.
Defaults to None.
full_support (tuple[int, ...] | None, optional): The qubits the returned tensor
will be defined over. Defaults to None for only using the qubit_support.
Returns:
Tensor: Output density matrix from solver.
"""
if not isinstance(state, DensityMatrix):
state = density_mat(state)
sol = mesolve(
H,
state,
self._noise_operators(full_support),
tsave,
solver,
options=options,
)
return DensityMatrix(sol.states[-1])

def __add__(self, other: AnalogNoise) -> AnalogNoise:
return AnalogNoise(
self.noise_operators + other.noise_operators,
self._qubit_support + other._qubit_support,
)


class AnalogDepolarizing(AnalogNoise):
"""
Defines jump operators for an Analog Depolarizing channel.
Under the depolarizing noise, a system in any state evolves
to the maximally mixed state at a rate :math:`p`.
The corresponding jump operators are :
.. math::
`L_{0,1,2} = \sqrt{\frac{p}{4}} \sigma_{x,y,z}`
where :math:`\sigma_{x,y,z}` correspond to the unitaries of the X,Y,Z gates.
Args:
error_param (float): Rate of depolarizing.
qubit_support (int | tuple[int, ...] | Support): Qubits defining the operation.
"""

def __init__(
self,
error_param: float,
qubit_support: int | tuple[int, ...] | Support,
) -> None:
"""Initializes AnalogDepolarizing.
Args:
error_param (float): Rate of depolarizing.
qubit_support (int | tuple[int, ...] | Support): Qubits defining the operation.
"""
coeff = sqrt(error_param / 4)
L0 = coeff * XMAT.squeeze()
L1 = coeff * YMAT.squeeze()
L2 = coeff * ZMAT.squeeze()
noise_operators: list[Tensor] = [L0, L1, L2]
super().__init__(noise_operators, qubit_support)
File renamed without changes.
Loading

0 comments on commit ef035ec

Please sign in to comment.