Hey there,
i am currently implementing a quantum variational algorithm, where my goal is to minimize a cost function, where i use a gradient based optimzer. Since i have no analytically closed form of my gradient, i am currently stuck.

I know one can use automatic differentiation when one has a pure state by
simply using the function expect’( op::AbstractBlock, reg).

Now my question: Is there a similar function for mixed states, where i have a density matrix, so i could simply use expect’(op::AbstractBlock, density_matrix)?

If that’s not the case: Do you have any idea how i could get the gradient of an expectation value in a simple fashion when my system is described by a density matrix in yao?

My expectation value of an Operator O is Tr[U+ O U rho] where the + stands for the hermitian conjugate and rho is the density matrix. The unitary time evolution operator U depends on a set of variational parameters. The gradient is calculated in respect to those parameters.

I think this feature is not available yet in Yao. Yao has a naive AD engine that only supports differentiating reversible gates, however density matrices simulation is usually not reversible and requires checkpointing.

Currently, I am focusing more on simulating quantum circuits by contracting tensor networks. The contraction is backended by OMEinsum, which has good support to automatic differentiation.

I am happy to help convert your simulation to tensor networks for better automatic differentiation if you could provide a minimum working example.

Thanks a lot for your answer and the offer. I really appreciate it. I will need to talk to my supervisor if that’s what we want.

One further question regarding automatic differentation, this time with a pure state.
Again there i am trying to get the gradient of my loss function in respect to my parameters.

Here a minimal working example:

using Yao
N = 2 #number of qubits
T = 0.1 #time step
Ψ = ArrayReg(complex([1.,0.,0.,0.])) #initial state
hamiltonian(x1,x2,x3) = x1*put(N, 1=>Z)*put(N, 2=>Z) + x2*put(N, 1=>X) + x3*put(N, 2=>X)
observable = put(N, 1=>Z)*put(N, 2=>Z)
C = expect(observable, Ψ => time_evolve(hamiltonian(1.,2.,3.), T)) #cost function

So what i want here is the gradient ∇ = (dC/dx1 , dC/dx2, dC/dx3)

but there is only one value saved in ∇. I assume its the derivative in respect to time.
How to i obtain the gradient of the loss function in respect to these parameters?

In your tasks, if errors are acceptable, a straight-forward approach is to trotterize the Hamiltonian and differentiate over the parameters with the built-in AD engine.

I can’t speak to Yao specifically, but fundamentally, there’s nothing different about open quantum systems vs. closed quantum systems. In closed quantum systems, you have |Ψ̇⟩ = H |Ψ⟩, in open quantum systems you have ρ̇ = L ρ. In both cases, Ψ/ρ are just vectors (vectorized density matrices: just write the columns of the density matrix underneath each other), and H/L are operators (super-operators, but again: if you’ve vectorized ρ, L is just a matrix). The only difference is that H has real eigenvalues (Hermitian) and L does not (or, only if there is no dissipation). Small detail: most textbooks put an extra imaginary-i in the Liouville equation, but I find it much more useful to have the exact same form of the Schrödinger and Liouville equation. I would point you to the docstring of liouvillian in QuantumControl.lj.

So if you have an analytic way to calculate the gradient for the closed case, you should also be able to calculate the gradient for the open case as well. Beyond that, I would also argue that auto-differentiating through a Hamiltonian/Liouvillian, respectively through a time propagation/circuit evaluation is a bad idea. I’ve proposed using “semi-automatic differentiation” which effectively avoids all of the AD overhead (assuming you can differentiate analytically the “standard” overlap functionals). I’ve implemented in this concept for pulse-level control. Yao operates on the higher circuit level, of course. I think the same idea applies in this context; but nobody has implemented semi-AD in Yao, I would assume. Still, you might find the paper instructive and be able to code something up yourself relatively easily.

Great, thanks for the hint. Trotterization and then using the AD-engine works fine. Since my gates share parameters, i needed to sum up the output of expect’.
Thanks for your help, appreciate your work a lot!