[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

@@ -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
 

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"]

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,

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,

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: