[Feature] Time-dependent Hamevo with noise operators (#621)

docs/content/time_dependent.md

@@ -1,4 +1,6 @@
-For use cases when the Hamiltonian of the system is time-dependent, Qadence provides a special parameter `TimePrameter("t")` that denotes the explicit time dependence. Using this time parameter one can define a parameterized block acting as the generator passed to `HamEvo` that encapsulates the required time dependence function.
+For use cases when the Hamiltonian of the system is time-dependent, Qadence provides a special parameter `TimeParameter("t")` that denotes the explicit time dependence. Using this time parameter, one can define a parameterized block acting as the generator passed to `HamEvo` that encapsulates the required time dependence function.
+
+# Noiseless time-dependent Hamiltonian evolution
 
 ```python exec="on" source="material-block" session="getting_started"
 from qadence import X, Y, HamEvo, TimeParameter, FeatureParameter, run
@@ -14,10 +16,10 @@ t = TimeParameter("t")
 omega_param = FeatureParameter("omega")
 
 # Arbitrarily compose a time-dependent generator
-generator_td = omega_param * (t * X(0) + t**2 * Y(1))
+td_generator = omega_param * (t * X(0) + t**2 * Y(1))
 
 # Create parameterized HamEvo block
-hamevo = HamEvo(generator_td, t)
+hamevo = HamEvo(td_generator, t)
 ```
 
 Note that when using `HamEvo` with a time-dependent generator, the actual time parameter that was used to construct the generator must be passed for the second argument `parameter`.
@@ -43,3 +45,39 @@ print(out_state)
 ```
 
 Note that Qadence makes no assumption on units. The unit of passed duration value $\tau$ must be aligned with the units of other parameters in the time-dependent generator so that the integral of generator $\overset{\tau}{\underset{0}{\int}}\mathcal{\hat{H}}(t){\rm d}t$ is dimensionless.
+
+# Noisy time-dependent Hamiltonian evolution
+
+To perform noisy time-dependent Hamiltonian evolution, one needs to pass a list of noise operators to the `noise_operators` argument in `HamEvo`. They correspond to the jump operators used within the time-dependent Schrodinger equation solver method `SolverType.DP5_ME`:
+
+```python exec="on" source="material-block" session="getting_started"
+from qadence import X, Y, HamEvo, TimeParameter, FeatureParameter, run
+from pyqtorch.utils import SolverType
+import torch
+
+# Simulation parameters
+ode_solver = SolverType.DP5_ME  # time-dependent Schrodinger equation solver method
+n_steps_hevo = 500  # integration time steps used by solver
+
+# Define block parameters
+t = TimeParameter("t")
+omega_param = FeatureParameter("omega")
+
+# Arbitrarily compose a time-dependent generator
+td_generator = omega_param * (t * X(0) + t**2 * Y(1))
+
+# Create parameterized HamEvo block
+noise_operators = [X(i) for i in td_generator.qubit_support]
+hamevo = HamEvo(td_generator, t, noise_operators = noise_operators)
+
+values = {"omega": torch.tensor(10.0), "duration": torch.tensor(1.0)}
+
+config = {"ode_solver": ode_solver, "n_steps_hevo": n_steps_hevo}
+
+out_state = run(hamevo, values = values, configuration = config)
+
+print(out_state)
+```
+
+!!! warning "Noise operators definition"
+    Note it is not possible to define `noise_operators` with parametric operators. If you want to do so, we recommend obtaining the tensors via run and set `noise_operators` using `MatrixBlock`. Also, `noise_operators` should have the same or a subset of the qubit support of the `HamEvo` instance.

mkdocs.yml

@@ -16,11 +16,11 @@ nav:
   - Contents:
     - Block system: content/block_system.md
     - Parametric programs: content/parameters.md
-    - Time-dependent generators: content/time_dependent.md
     - Quantum models: content/quantummodels.md
     - Quantum registers: content/register.md
     - State initialization: content/state_init.md
     - Arbitrary Hamiltonians: content/hamiltonians.md
+    - Time-dependent generators: content/time_dependent.md
     - QML Constructors: content/qml_constructors.md
     - Wavefunction overlaps: content/overlap.md
     - Backends: content/backends.md

pyproject.toml

@@ -26,7 +26,7 @@ authors = [
 ]
 requires-python = ">=3.9"
 license = { text = "Apache 2.0" }
-version = "1.8.0"
+version = "1.9.0"
 classifiers = [
     "License :: OSI Approved :: Apache Software License",
     "Programming Language :: Python",

qadence/backends/pyqtorch/convert_ops.py

@@ -245,13 +245,26 @@ def convert_block(
             generator = convert_block(block.generator, n_qubits, config)[0]  # type: ignore[arg-type]
         time_param = config.get_param_name(block)[0]
 
+        # convert noise operators here
+        noise_operators: list = [
+            convert_block(noise_block, config=config)[0] for noise_block in block.noise_operators
+        ]
+        if len(noise_operators) > 0:
+            # squeeze batch size for noise operators
+            noise_operators = [
+                pyq_op.tensor(full_support=qubit_support).squeeze(-1) for pyq_op in noise_operators
+            ]
+
         return [
             pyq.HamiltonianEvolution(
                 qubit_support=qubit_support,
                 generator=generator,
                 time=time_param,
                 cache_length=0,
                 duration=duration,
+                solver=config.ode_solver,
+                steps=config.n_steps_hevo,
+                noise_operators=noise_operators,
             )
         ]
 

qadence/blocks/primitive.py

@@ -356,6 +356,7 @@ class TimeEvolutionBlock(ParametricBlock):
     """
 
     name = "TimeEvolutionBlock"
+    noise_operators: list = list()
 
     @property
     def has_parametrized_generator(self) -> bool:

qadence/operations/ham_evo.py

@@ -55,6 +55,8 @@ class HamEvo(TimeEvolutionBlock):
         duration: (optional) duration of the evolution in case of time-dependent
             generator. By default, a FeatureParameter with tag "duration" will
             be initialized, and the value will then be required in the values dict.
+        noise_operators: (optional) the list of jump operators to use when using
+            a shrodinger solver, allowing to perform noisy simulations.
 
     Examples:
 
@@ -80,6 +82,10 @@ class HamEvo(TimeEvolutionBlock):
     hamiltonian = t * add(X(i) for i in range(n_qubits))
     hevo = HamEvo(hamiltonian, parameter=t)
     state = run(hevo, values = {"duration": torch.tensor(1.0)})
+
+    # Adding noise operators
+    noise_ops = [X(0)]
+    hevo = HamEvo(hamiltonian, parameter=t, noise_operators=noise_ops)
     ```
     """
 
@@ -92,6 +98,7 @@ def __init__(
         parameter: TParameter,
         qubit_support: tuple[int, ...] = None,
         duration: TParameter | None = None,
+        noise_operators: list[AbstractBlock] = list(),
     ):
         params = {}
         if qubit_support is None and not isinstance(generator, AbstractBlock):
@@ -147,6 +154,21 @@ def __init__(
         self.generator = generator
         self.duration = duration
 
+        if len(noise_operators) > 0:
+            if not all(
+                [
+                    len(set(op.qubit_support + self.qubit_support) - set(self.qubit_support)) == 0
+                    for op in noise_operators
+                ]
+            ):
+                raise ValueError(
+                    "Noise operators should be defined"
+                    " over the same or a subset of the qubit support"
+                )
+            if True in [op.is_parametric for op in noise_operators]:
+                raise ValueError("Parametric operators are not supported")
+        self.noise_operators = noise_operators
+
     @classmethod
     def num_parameters(cls) -> int:
         return 2

tests/backends/pyq/test_time_dependent_generator.py

@@ -9,10 +9,20 @@
 from metrics import MIDDLE_ACCEPTANCE
 from pyqtorch.utils import SolverType
 
-from qadence import AbstractBlock, HamEvo, QuantumCircuit, QuantumModel, Register, run
+from qadence import (
+    AbstractBlock,
+    HamEvo,
+    QuantumCircuit,
+    QuantumModel,
+    Register,
+    block_to_tensor,
+    run,
+)
+from qadence.blocks import MatrixBlock
+from qadence.operations import RZ, I, X
 
 
-@pytest.mark.parametrize("duration", [0.5, 1.0, 2.0])
+@pytest.mark.parametrize("duration", [0.5, 1.0])
 @pytest.mark.parametrize("ode_solver", [SolverType.DP5_SE, SolverType.KRYLOV_SE])
 def test_time_dependent_generator(
     qadence_generator: AbstractBlock,
@@ -25,12 +35,18 @@ def test_time_dependent_generator(
 ) -> None:
     n_steps = 500
 
+    psi_0 = qutip.basis(4, 0)
+
     # simulate with qadence HamEvo using QuantumModel
     hamevo = HamEvo(qadence_generator, time_param)
     reg = Register(2)
     circ = QuantumCircuit(reg, hamevo)
     config = {"ode_solver": ode_solver, "n_steps_hevo": n_steps}
-    values = values = {"x": torch.tensor(feature_param_x), "duration": torch.tensor(duration)}
+    values = {
+        "x": torch.tensor(feature_param_x),
+        "y": torch.tensor(feature_param_y),
+        "duration": torch.tensor(duration),
+    }
 
     model = QuantumModel(circ, configuration=config)
     state_qadence0 = model.run(values=values)
@@ -40,9 +56,82 @@ def test_time_dependent_generator(
 
     # simulate with qutip
     t_points = np.linspace(0, duration, n_steps)
-    psi_0 = qutip.basis(4, 0)
+
     result = qutip.sesolve(qutip_generator, psi_0, t_points)
+
     state_qutip = torch.tensor(result.states[-1].full().T)
 
     assert torch.allclose(state_qadence0, state_qutip, atol=MIDDLE_ACCEPTANCE)
     assert torch.allclose(state_qadence1, state_qutip, atol=MIDDLE_ACCEPTANCE)
+
+
+@pytest.mark.parametrize("duration", [0.5, 1.0])
+@pytest.mark.parametrize(
+    "noise_op",
+    [
+        I(0) * I(1),
+        X(0),
+        MatrixBlock(
+            block_to_tensor(X(0), qubit_support=tuple(range(2)), use_full_support=True),
+            qubit_support=tuple(range(2)),
+        ),
+    ],
+)
+def test_noisy_time_dependent_generator(
+    qadence_generator: AbstractBlock,
+    qutip_generator: Callable,
+    time_param: str,
+    feature_param_x: float,
+    feature_param_y: float,
+    duration: float,
+    noise_op: AbstractBlock,
+) -> None:
+    n_steps = 500
+    ode_solver = SolverType.DP5_ME
+    n_qubits = 2
+
+    # Define jump operators
+    list_ops = [noise_op]
+
+    # simulate with qadence HamEvo using QuantumModel
+    hamevo = HamEvo(qadence_generator, time_param, noise_operators=list_ops)
+    reg = Register(n_qubits)
+    circ = QuantumCircuit(reg, hamevo)
+    n_qubits = circ.n_qubits
+
+    config = {"ode_solver": ode_solver, "n_steps_hevo": n_steps}
+    values = {
+        "x": torch.tensor(feature_param_x),
+        "y": torch.tensor(feature_param_y),
+        "duration": torch.tensor(duration),
+    }
+
+    model = QuantumModel(circ, configuration=config)
+    state_qadence0 = model.run(values=values)
+
+    # simulate with qadence.execution
+    state_qadence1 = run(hamevo, values=values, configuration=config)
+
+    # simulate with qutip
+    t_points = np.linspace(0, duration, n_steps)
+    noise_tensor = (
+        block_to_tensor(noise_op, qubit_support=tuple(range(n_qubits)), use_full_support=True)
+        .squeeze(0)
+        .numpy()
+    )
+    list_ops_qutip = [qutip.Qobj(noise_tensor)]
+    result = qutip.mesolve(qutip_generator, qutip.basis(2**n_qubits, 0), t_points, list_ops_qutip)
+
+    state_qutip = torch.tensor(result.states[-1].full()).unsqueeze(0)
+    assert torch.allclose(state_qadence0, state_qutip, atol=MIDDLE_ACCEPTANCE)
+    assert torch.allclose(state_qadence1, state_qutip, atol=MIDDLE_ACCEPTANCE)
+
+
+@pytest.mark.parametrize("noise_op", [I(0) * I(1) * I(3), X(3), RZ(0, "theta")])
+def test_error_noise_operators_hamevo(
+    qadence_generator: AbstractBlock,
+    time_param: str,
+    noise_op: AbstractBlock,
+) -> None:
+    with pytest.raises(ValueError):
+        hamevo = HamEvo(qadence_generator, time_param, noise_operators=[noise_op])

tests/conftest.py

@@ -268,9 +268,12 @@ def qadence_generator(omega: float, time_param: str, feature_param_y: float) ->
 @fixture
 def qutip_generator(omega: float, feature_param_x: float, feature_param_y: float) -> Callable:
     def generator_t(t: float, args: Any) -> qutip.Qobj:
-        return omega * (
-            np.sin(feature_param_y * t) * qutip.tensor(qutip.sigmax(), qutip.qeye(2))
-            + feature_param_x * t**2 * qutip.tensor(qutip.qeye(2), qutip.sigmay())
+        return qutip.Qobj(
+            omega
+            * (
+                np.sin(feature_param_y * t) * qutip.tensor(qutip.sigmax(), qutip.qeye(2))
+                + feature_param_x * (t**2) * qutip.tensor(qutip.qeye(2), qutip.sigmay())
+            ).full()
         )
 
     return generator_t