[unitaryhack] Add closed- and open-system dynamics propagators

[friedsam]

Closes #3437.

Adds cudaq.contrib.propagator for closed- and open-system dynamics.

For closed-system dynamics, the helper computes the propagator U(t) satisfying
dU/dt = -i H(t) U, U(t_initial) = I
by evolving the identity matrix under the left-multiplication superoperator generated by -i H(t).

For open-system dynamics with collapse operators, the helper computes the Lindblad superoperator propagator S(t) satisfying
d vec(rho)/dt = L(t) vec(rho), S(t_initial) = I
by evolving the Liouville basis states through cudaq.evolve with the supplied collapse operators and stacking the evolved basis states into the dense propagator matrix.

Implementation:

• adds cudaq.contrib.propagator
• supports closed-system propagators via cudaq.SuperOperator.left_multiply(-1j * H)
• supports open-system Lindblad superoperator propagators via collapse operators passed to cudaq.evolve
• supports final propagator output
• supports intermediate propagators along the schedule
• supports batched Hamiltonians through cudaq.evolve
• supports shared and per-Hamiltonian collapse-operator lists for batched open-system runs
• returns dense NumPy matrices:
  • (dim, dim) for closed-system propagators
  • (dim**2, dim**2) for open-system Lindblad propagators

Tests:

• constant Z Hamiltonian vs scipy.linalg.expm
• constant X Hamiltonian vs scipy.linalg.expm
• batched Hamiltonians
• intermediate closed-system propagators
• commuting time-dependent Hamiltonian with analytic integral
• noncommuting time-dependent Hamiltonian vs scipy.integrate.solve_ivp
• multi-degree-of-freedom propagator shape
• open-system amplitude damping propagator vs scipy.linalg.expm of the Lindblad generator
• open-system propagator shape
• open-system schedule parameter name handling
• open-system nonzero Hamiltonian
• batched open-system propagators with shared collapse operators
• batched open-system propagators with per-Hamiltonian collapse-operator lists
• intermediate open-system propagators

Local validation:

• python -m yapf -i python/cudaq/contrib/propagator.py
• python -m py_compile python/cudaq/contrib/propagator.py python/tests/contrib/test_propagator.py
• python -m pytest python/tests/contrib/test_propagator.py -q -rs
• python -m pre_commit run spellcheck --all-files --hook-stage pre-push
• python -m pre_commit run --files python/cudaq/contrib/propagator.py
• git diff --check

GPU validation:

• Brev L40S
• python -m pytest python/tests/contrib/test_propagator.py -q -vv -rs
• 14 passed in 2.38s

AI assistance disclosure:
I used AI to help reason through the CUDA-Q dynamics API and to help with drafting/debugging for the initial implementation and tests. The numerical tests are GPU-gated and were skipped locally because no NVIDIA GPU/dynamics backend is available on my machine.

[copy-pr-bot]

This pull request requires additional validation before any workflows can run on NVIDIA's runners.

Pull request vetters can view their responsibilities here.

Contributors can view more details about this message here.

[1tnguyen]

Thanks @friedsam for the contribution!

Could you also please extend this PR to cover open-quantum systems (with collapse operators) as mentioned in #3437 (comment)?

For example, the propagator ($S$) of the Lindblad master equation:

$$ \frac{d\rho}{dt} = -i[H(t),\rho] + \sum_k \left( C_k(t)\rho C_k^\dagger(t) -\frac{1}{2}C_k^\dagger(t)C_k(t)\rho -\frac{1}{2}\rho C_k^\dagger(t)C_k(t) \right) $$

would follow this ODE:

$$ \frac{dS}{dt} = \left[ -i\left(I_d \otimes H\right) +i\left(H^T \otimes I_d\right) + \sum_k \left( C_k^* \otimes C_k -\frac{1}{2} I_d \otimes C_k^\dagger C_k -\frac{1}{2} \left(C_k^\dagger C_k\right)^T \otimes I_d \right) \right]S \qquad $$

[friedsam]

Thanks, yes. I’ll extend the helper to support collapse operators and return the Lindblad superoperator propagator for open-system dynamics.

[1tnguyen]

/ok to test f527c3d

Command Bot: Processing...

[friedsam]

Thanks. I extended the helper to support collapse operators and return the Lindblad superoperator propagator for open-system dynamics. The implementation follows the existing CUDA-Q SuperOperator construction pattern using left_multiply, right_multiply, and left_right_multiply.

Local validation:

• python -m py_compile python/cudaq/contrib/propagator.py python/tests/contrib/test_propagator.py
• python -m pytest python/tests/contrib/test_propagator.py -q -rs
• git diff --check

[1tnguyen]

Hi @friedsam,

The test is failing:

============================= test session starts ==============================
  platform linux -- Python 3.12.3, pytest-9.0.3, pluggy-1.6.0
  rootdir: /tmp/cuda-quantum-pr4672
  configfile: pyproject.toml
  plugins: anyio-4.13.0
  collected 8 items

  python/tests/contrib/test_propagator.py .....FFF     [100%]

For this task, I think you'd need a machine with NVIDIA GPU to do the validation/debug.

[friedsam]

Thanks. Could you please paste the traceback for the three failing tests? I can likely patch from the failure details. I’m also looking into a temporary NVIDIA GPU environment, but the traceback would be the fastest path.

[1tnguyen]

Thanks. Could you please paste the traceback for the three failing tests? I can likely patch from the failure details. I’m also looking into a temporary NVIDIA GPU environment, but the traceback would be the fastest path.

Here is the log:

 ============================= test session starts ==============================
  platform linux -- Python 3.12.3, pytest-9.0.3, pluggy-1.6.0
  rootdir: /tmp/cuda-quantum-pr4672
  configfile: pyproject.toml
  plugins: anyio-4.13.0
  collected 8 items

  cuda-quantum-pr4672/python/tests/contrib/test_propagator.py .....FFF     [100%]

  =================================== FAILURES ===================================
  ___________ test_propagator_noncommuting_time_dependent_hamiltonian ____________

      def test_propagator_noncommuting_time_dependent_hamiltonian():
          ax = 0.2
          bz = 0.35
          t_final = 0.8

          hamiltonian = (ax * ScalarOperator(lambda t: t) * spin.x(0) +
                         bz * ScalarOperator(lambda t: 1.0 - t) * spin.z(0))
          schedule = Schedule(np.linspace(0.0, t_final, 41), ["t"])

          computed = cudaq.contrib.propagator(hamiltonian, {0: 2}, schedule)

          def hamiltonian_matrix(t):
              return ax * t * _X + bz * (1.0 - t) * _Z

          expected = _solve_propagator_reference(hamiltonian_matrix, t_final)

  >       np.testing.assert_allclose(computed, expected, atol=1e-4)
  E       AssertionError:
  E       Not equal to tolerance rtol=1e-07, atol=0.0001
  E
  E       Mismatched elements: 2 / 4 (50%)
  E       Max absolute difference among violations: 0.01190243
  E       Max relative difference among violations: 0.1865569
  E        ACTUAL: array([[ 0.983887-0.16702j , -0.005951-0.063522j],
  E              [ 0.005951-0.063522j,  0.983887+0.16702j ]])
  E        DESIRED: array([[ 0.983887-0.16702j ,  0.005951-0.063522j],
  E              [-0.005951-0.063522j,  0.983887+0.16702j ]])

  cuda-quantum-pr4672/python/tests/contrib/test_propagator.py:152: AssertionError
  ________________ test_open_system_propagator_amplitude_damping _________________

      def test_open_system_propagator_amplitude_damping():
          gamma = 0.2
          t_final = 0.7

          hamiltonian = 0.0 * spin.z(0)
          collapse_operators = [np.sqrt(gamma) * spin.minus(0)]
          schedule = Schedule(np.linspace(0.0, t_final, 21), ["t"])

  >       computed = cudaq.contrib.propagator(
              hamiltonian,
              {0: 2},
              schedule,
              collapse_operators=collapse_operators,
          )

  cuda-quantum-pr4672/python/tests/contrib/test_propagator.py:177:
  _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
  cudaq-pr4672-install/cudaq/contrib/propagator.py:113: in propagator
      _open_system_generator(h, collapse_operators) for h in hamiltonians
  _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

  hamiltonian = <cudaq.mlir._mlir_libs._quakeDialects.cudaq_runtime.SpinOperatorTerm object at 0x7548451a7c60>
  collapse_operators = [<cudaq.mlir._mlir_libs._quakeDialects.cudaq_runtime.SpinOperator object at 0x7548451a7d80>]

      def _open_system_generator(hamiltonian, collapse_operators):
          import cudaq

          generator = cudaq.SuperOperator()
          generator += cudaq.SuperOperator.left_multiply(-1j * hamiltonian)
          generator += cudaq.SuperOperator.right_multiply(1j * hamiltonian)

          for collapse_operator in collapse_operators:
  >           collapse_operator_dagger = collapse_operator.dagger()
  E           AttributeError: 'cudaq.mlir._mlir_libs._quakeDialects.cudaq_runtime.SpinOperator' object has no attribute 'dagger'

  cudaq-pr4672-install/cudaq/contrib/propagator.py:51: AttributeError
  ______________________ test_open_system_propagator_shape _______________________

      def test_open_system_propagator_shape():
          gamma = 0.1
          t_final = 0.3

          hamiltonian = 0.0 * spin.z(0)
          collapse_operators = [np.sqrt(gamma) * spin.minus(0)]
          schedule = Schedule(np.linspace(0.0, t_final, 11), ["t"])

  >       computed = cudaq.contrib.propagator(
              hamiltonian,
              {0: 2},
              schedule,
              collapse_operators=collapse_operators,
          )

  cuda-quantum-pr4672/python/tests/contrib/test_propagator.py:198:
  _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
  cudaq-pr4672-install/cudaq/contrib/propagator.py:113: in propagator
      _open_system_generator(h, collapse_operators) for h in hamiltonians
  _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

  hamiltonian = <cudaq.mlir._mlir_libs._quakeDialects.cudaq_runtime.SpinOperatorTerm object at 0x7548451a7e10>
  collapse_operators = [<cudaq.mlir._mlir_libs._quakeDialects.cudaq_runtime.SpinOperator object at 0x7548451a7f30>]

      def _open_system_generator(hamiltonian, collapse_operators):
          import cudaq

          generator = cudaq.SuperOperator()
          generator += cudaq.SuperOperator.left_multiply(-1j * hamiltonian)
          generator += cudaq.SuperOperator.right_multiply(1j * hamiltonian)

          for collapse_operator in collapse_operators:
  >           collapse_operator_dagger = collapse_operator.dagger()
  E           AttributeError: 'cudaq.mlir._mlir_libs._quakeDialects.cudaq_runtime.SpinOperator' object has no attribute 'dagger'

  cudaq-pr4672-install/cudaq/contrib/propagator.py:51: AttributeError
  =========================== short test summary info ============================
  FAILED cuda-quantum-pr4672/python/tests/contrib/test_propagator.py::test_propagator_noncommuting_time_dependent_hamiltonian
  FAILED cuda-quantum-pr4672/python/tests/contrib/test_propagator.py::test_open_system_propagator_amplitude_damping
  FAILED cuda-quantum-pr4672/python/tests/contrib/test_propagator.py::test_open_system_propagator_shape
  ========================= 3 failed, 5 passed in 1.54s ==========================

[github-actions]

CI Summary (push) — ✅ passed

Run #27253771876 · ✅ 6 · ⏩ 7 · ❌ 0 · ⛔ 0

Top-level jobs (13)

Job	Result
binaries	⏩ skipped
build_and_test	✅ success
config_devdeps	✅ success
config_source_build	⏩ skipped
config_wheeldeps	✅ success
devdeps	✅ success
docker_image	⏩ skipped
gen_code_coverage	⏩ skipped
metadata	✅ success
python_metapackages	⏩ skipped
python_wheels	⏩ skipped
source_build	⏩ skipped
wheeldeps	✅ success

⏩ Skipped jobs (7) — intentionally skipped on PR builds; run on merge_group / workflow_dispatch

Job
binaries
config_source_build
docker_image
gen_code_coverage
python_metapackages
python_wheels
source_build

All sub-jobs (42) — every matrix leg, with links

Job	Status	Link
Build and test (amd64, gcc12, openmpi) / Dev environment (Debug)	✅ success	view
Build and test (amd64, gcc12, openmpi) / Dev environment (Python)	✅ success	view
Build and test (amd64, llvm, openmpi) / Dev environment (Debug)	✅ success	view
Build and test (amd64, llvm, openmpi) / Dev environment (Python)	✅ success	view
Build and test (arm64, llvm, openmpi) / Dev environment (Debug)	✅ success	view
Build and test (arm64, llvm, openmpi) / Dev environment (Python)	✅ success	view
CI Summary	❔ in_progress	view
Configure build (devdeps)	✅ success	view
Configure build (source_build)	⏩ skipped	view
Configure build (wheeldeps)	✅ success	view
Create CUDA Quantum installer	⏩ skipped	view
Create Docker images	⏩ skipped	view
Create Python metapackages	⏩ skipped	view
Create Python wheels	⏩ skipped	view
Gen code coverage	⏩ skipped	view
Load dependencies (amd64, gcc12) / Caching	✅ success	view
Load dependencies (amd64, gcc12) / Finalize	✅ success	view
Load dependencies (amd64, gcc12) / Metadata	✅ success	view
Load dependencies (amd64, llvm) / Caching	✅ success	view
Load dependencies (amd64, llvm) / Finalize	✅ success	view
Load dependencies (amd64, llvm) / Metadata	✅ success	view
Load dependencies (arm64, gcc12) / Caching	✅ success	view
Load dependencies (arm64, gcc12) / Finalize	✅ success	view
Load dependencies (arm64, gcc12) / Metadata	✅ success	view
Load dependencies (arm64, llvm) / Caching	✅ success	view
Load dependencies (arm64, llvm) / Finalize	✅ success	view
Load dependencies (arm64, llvm) / Metadata	✅ success	view
Load source build cache	⏩ skipped	view
Load wheel dependencies (amd64, 12.6) / Caching	✅ success	view
Load wheel dependencies (amd64, 12.6) / Finalize	✅ success	view
Load wheel dependencies (amd64, 12.6) / Metadata	✅ success	view
Load wheel dependencies (amd64, 13.0) / Caching	✅ success	view
Load wheel dependencies (amd64, 13.0) / Finalize	✅ success	view
Load wheel dependencies (amd64, 13.0) / Metadata	✅ success	view
Load wheel dependencies (arm64, 12.6) / Caching	✅ success	view
Load wheel dependencies (arm64, 12.6) / Finalize	✅ success	view
Load wheel dependencies (arm64, 12.6) / Metadata	✅ success	view
Load wheel dependencies (arm64, 13.0) / Caching	✅ success	view
Load wheel dependencies (arm64, 13.0) / Finalize	✅ success	view
Load wheel dependencies (arm64, 13.0) / Metadata	✅ success	view
Prepare cache clean-up	❔ in_progress	view
Retrieve PR info	✅ success	view

✅ Required checks (6/6) — declared in .github/required-checks.yml for push

Required check	Status	Link
Build and test (amd64, llvm, openmpi) / Dev environment (Debug)	✅ success	view
Build and test (amd64, llvm, openmpi) / Dev environment (Python)	✅ success	view
Build and test (arm64, llvm, openmpi) / Dev environment (Debug)	✅ success	view
Build and test (arm64, llvm, openmpi) / Dev environment (Python)	✅ success	view
Build and test (amd64, gcc12, openmpi) / Dev environment (Debug)	✅ success	view
Build and test (amd64, gcc12, openmpi) / Dev environment (Python)	✅ success	view

[friedsam]

Thanks, that makes sense. I see the distinction: the current open-system branch constructs the Lindblad RHS for density-matrix evolution, while the requested object is the propagator S satisfying dS/dt = L(t) S.

I’ll revise the open-system implementation to construct the Liouvillian on vec(rho) and evolve the identity in that enlarged space.

[friedsam]

Thanks — I think the API-compatible way to implement the open-system propagator is to assemble S column-by-column:

• build the Lindblad SuperOperator L for d vec(rho)/dt = L vec(rho)
• evolve each Liouville basis vector e_j under L
• stack the evolved basis vectors as columns of S

This should give the requested map vec(rho(t)) = S vec(rho(0)), equivalent to solving dS/dt = L(t) S, without requiring a separate dense Liouvillian operator constructor.

Does that match what you had in mind?

[friedsam]

Updated. The open-system branch now assembles the propagator map by evolving the Liouville-space basis vectors under the Lindblad SuperOperator and stacking the evolved vectors as columns of S, so the returned matrix satisfies vec(rho(t)) = S vec(rho(0)).

Local checks pass. The pytest module remains GPU-gated and skips locally.

[friedsam]

You are very correct with this. I’ll switch this from a Python-side serial loop to a batched cudaq.evolve call over the Liouville basis states, so the columns are evolved as a batch rather than one-by-one.

[1tnguyen]

Thanks — I think the API-compatible way to implement the open-system propagator is to assemble S column-by-column:

• build the Lindblad SuperOperator L for d vec(rho)/dt = L vec(rho)
• evolve each Liouville basis vector e_j under L
• stack the evolved basis vectors as columns of S

This should give the requested map vec(rho(t)) = S vec(rho(0)), equivalent to solving dS/dt = L(t) S, without requiring a separate dense Liouvillian operator constructor.

Does that match what you had in mind?

I posted a comment in the code. I think this is not an efficient way to use the GPU as we made the computation sequential.
You can take a look at the comment, which has a reference for the propagator ODE.

For testing coverage, we can test normal density matrix state evolve and a propagator calculation based on the same dynamics. The final density matrix $\rho$ should match the direct application of the propagator on the initial state.

[1tnguyen]

Local checks pass. The pytest module remains GPU-gated and skips locally.

Thanks! Please let us know when you get a chance to actually validate all the tests so that we can take a closer look at the implementation.

[friedsam]

Updated and GPU-tested.

I changed the open-system implementation to construct the dense Liouvillian L(t) directly and evolve the identity in Liouville space with dS/dt = L(t) S, instead of assembling S by evolving basis columns.

I also ran the contrib propagator tests on a Brev NVIDIA L40S GPU instance:

===================================================== test session starts ======================================================
platform linux -- Python 3.10.12, pytest-9.0.3, pluggy-1.6.0 -- /home/ubuntu/cuda-quantum/cudaq-test/bin/python
cachedir: .pytest_cache
rootdir: /home/ubuntu/cuda-quantum
configfile: pyproject.toml
collected 8 items

python/tests/contrib/test_propagator.py::test_propagator_constant_z_hamiltonian PASSED [ 12%]
python/tests/contrib/test_propagator.py::test_propagator_constant_x_hamiltonian PASSED [ 25%]
python/tests/contrib/test_propagator.py::test_propagator_batched_hamiltonians PASSED [ 37%]
python/tests/contrib/test_propagator.py::test_propagator_intermediate_results PASSED [ 50%]
python/tests/contrib/test_propagator.py::test_propagator_commuting_time_dependent_hamiltonian PASSED [ 62%]
python/tests/contrib/test_propagator.py::test_propagator_noncommuting_time_dependent_hamiltonian PASSED [ 75%]
python/tests/contrib/test_propagator.py::test_open_system_propagator_amplitude_damping PASSED [ 87%]
python/tests/contrib/test_propagator.py::test_open_system_propagator_shape PASSED [100%]

====================================================== 8 passed in 1.42s =======================================================

[friedsam]

Updated to preserve the structured SuperOperator representation for open-system propagators instead of expanding collapse operators into dense full-system matrices.

This now builds the Lindblad generator using left_multiply, right_multiply, and left_right_multiply, and reconstructs the propagator by evolving Liouville basis states as a batch.

GPU-tested on Brev NVIDIA L40S:

python -m pytest python/tests/contrib/test_propagator.py -q -vv -rs

8 passed in 1.45s

[1tnguyen]

/ok to test a1c0ecf

Command Bot: Processing...

[1tnguyen]

Hi @friedsam,

Just checking in on this PR. Thanks again for opening it; this is a nice contribution, and it looks very close.

The remaining feedback is minor, mostly around code comments/test coverage.
If you have time to address those, I think we can get this ready to merge. And if anything is confusing or you’re unsure how to make the changes, please ask here. We're happy to help.

Really appreciate your contribution to CUDA-Q!

[friedsam]

Hi, thank you so much for all the feedback. I actually worked on it last evening but didn't push as I didn't want to disrupt the tests. I will finish it up tonight, test on brev again and commit everything

[friedsam]

Updated the batched open-system path to flatten the Hamiltonian/basis-state combinations before calling cudaq.evolve, then regroup the returned basis evolutions into one propagator per Hamiltonian. This avoids passing nested initial-state lists to the dynamics backend.

I also added the requested helper comments and tests for nonzero Hamiltonian plus collapse operators, shared collapse operators across Hamiltonian batches, and per-Hamiltonian collapse-operator lists.

GPU-tested on Brev NVIDIA L40S:

(cudaq-test) ubuntu@crusoe-prod--inst-3evbehfhoujpidcjh56yh7x5y5i:~/cuda-quantum$ python -m py_compile python/cudaq/contrib/propagator.py python/tests/contrib/test_propagator.py
python -m py_compile cudaq-test/lib/python3.10/site-packages/cudaq/contrib/propagator.py
python -m pytest python/tests/contrib/test_propagator.py -q -vv -rs
===================================================== test session starts ======================================================
platform linux -- Python 3.10.12, pytest-9.0.3, pluggy-1.6.0 -- /home/ubuntu/cuda-quantum/cudaq-test/bin/python
cachedir: .pytest_cache
rootdir: /home/ubuntu/cuda-quantum
configfile: pyproject.toml
collected 13 items

python/tests/contrib/test_propagator.py::test_propagator_constant_z_hamiltonian PASSED [ 7%]
python/tests/contrib/test_propagator.py::test_propagator_constant_x_hamiltonian PASSED [ 15%]
python/tests/contrib/test_propagator.py::test_propagator_batched_hamiltonians PASSED [ 23%]
python/tests/contrib/test_propagator.py::test_propagator_intermediate_results PASSED [ 30%]
python/tests/contrib/test_propagator.py::test_propagator_commuting_time_dependent_hamiltonian PASSED [ 38%]
python/tests/contrib/test_propagator.py::test_propagator_noncommuting_time_dependent_hamiltonian PASSED [ 46%]
python/tests/contrib/test_propagator.py::test_open_system_propagator_amplitude_damping PASSED [ 53%]
python/tests/contrib/test_propagator.py::test_open_system_propagator_shape PASSED [ 61%]
python/tests/contrib/test_propagator.py::test_propagator_multi_degree_of_freedom PASSED [ 69%]
python/tests/contrib/test_propagator.py::test_open_system_propagator_schedule_parameter_name PASSED [ 76%]
python/tests/contrib/test_propagator.py::test_open_system_propagator_nonzero_hamiltonian PASSED [ 84%]
python/tests/contrib/test_propagator.py::test_open_system_propagator_batched_shared_collapse_operators PASSED [ 92%]
python/tests/contrib/test_propagator.py::test_open_system_propagator_batched_collapse_operator_lists PASSED [100%]

====================================================== 13 passed in 1.75s ======================================================

[1tnguyen]

/ok to test e3ee81c

Command Bot: Processing...

[1tnguyen]

Re: https://github.com/NVIDIA/cuda-quantum/actions/runs/27244073636/job/80455376958?pr=4672

We can use backsticks (`) to escape those specialized terms in the Python code comments.

[1tnguyen]

/ok to test ba05bf2

Command Bot: Processing...

[1tnguyen]

@friedsam

BTW, if this PR is for unitaryHACK, please put [unitaryhack] in the title of the PR.

Also, could you please post a comment to the issue so that I can assign it to you?

Thanks!

[1tnguyen]

/ok to test 7fa78d6

Command Bot: Processing...

[1tnguyen]

/ok to test 1e43827

Command Bot: Processing...