Fidelity reported by GRAPE doesn't match manually computed one

[pmenczel]

The following script sets up a simple example problem where a (dissipative) qubit should be transferred from the ground state to the excited state. We use GRAPE and ask for an infidelity smaller then 0.001. We then use the optimized control returned by GRAPE to re-run the time evolution manually.

The final infidelity claimed by GRAPE is 0.0008. The infidelity found in the manual run is 0.0038.

A small discrepancy would be fine due to numerical inaccuracies, different integration methods etc, and small discrepancies also occur using other algorithms. However, when we ask for an infidelity of 0.001 and the result returned by GRAPE misses it by almost 300%, it seems like an issue to me.

import numpy as np
import qutip as qt
from qutip_qoc import Objective, optimize_pulses

def infidelity(state, target_state):
    """
    Infidelity used for density matrices in qutip-qtrl and qutip-qoc
    """
    diff = state - target_state
    return np.real(diff.overlap(diff)) / (2 * target_state.norm())

# --- Setup problem ---
Hd = qt.Qobj(np.diag([1, 2]))
c_ops = [np.sqrt(0.1) * qt.sigmam()]
Hc = qt.sigmax()

Ld = qt.liouvillian(H=Hd, c_ops=c_ops)
Lc = qt.liouvillian(Hc)
L = [Ld, Lc]

initial_state = qt.fock_dm(2, 0)
target_state = qt.fock_dm(2, 1)

times = np.linspace(0, 2 * np.pi, 250)

# --- Run GRAPE ---
guess_pulse = np.sin(times)

res_grape = optimize_pulses(
    objectives = Objective(initial_state, L, target_state),
    control_parameters = {
        "ctrl_x": {"guess": guess_pulse, "bounds": [-1, 1]},
    },
    tlist = times,
    algorithm_kwargs = {
        "alg": "GRAPE",
        "fid_err_targ": 0.001,
    },
)

grape_infidelity = res_grape.infidelity

# --- Manually run evolution ---
state = initial_state
dt = times[1] - times[0]

for c in res_grape.optimized_controls[0][:-1]:
    propagator = ((Ld + c * Lc) * dt).expm()
    state = propagator(state)

real_infidelity = infidelity(state, target_state)

# --- These should be the same ---
print(grape_infidelity)
print(real_infidelity)
assert np.isclose(grape_infidelity, real_infidelity, atol=1e-3)

Note that I do the manual evolution by just exponentiating the propagator on every time step, since GRAPE protocols are piecewise constant. Doing the manual evolution with mesolve gives the same result as this.

This issue was discovered while creating the interactive test notebooks in #44

[ajgpitch]

The relationship between the density matrix fidelity and the infidelity used within Qtrl is explained in this very handy document that I found.

[rochisha0]

As @nwlambert suggested, that qutip-qtrl adds a scale_factor to the fidelity. So now what would be the best way to fix this?

Should we manually compute the infidelity at the end for grape? @ajgpitch

[nwlambert]

i think there is still some issue in addition to this factor. strangely also the final state in res_grape.final_states[0] gives a slightly different result from grape_infidelity and real_infidelity.

I think probably this is all due to small differences in how the integration is done, and some resulting small errors/deviations. For example, the final_state object is calculated in qoc here (I think):

qutip-qoc/src/qutip_qoc/result.py

Line 327 in d5ff028

def final_states(self):

but this takes the time-dependent optimized hamiltonian, and finds the final_state by calling mesolve with tlist = [0,final_time]. This probably has intrinsically some numerical error due to limited time-steps (If i copy paste that code and put in more the 250 time-stepped tlist that we gave optimize_pulses, i get something a bit different again from everything else!).

meanwhile, the (H*dt).expm() in the qutip code above is probably a bit closer to what qtrl itself does internally, but not exactly equivalent.

so I guess, two things need to be done

• for state transfer, correct scale_factor in qtrl or the qoc interface (probably easier in the latter, since rochisha already added state_transfer checks? will qtrl-only users care about this?)

• for final_state, can we use the actual result returned by qtrl for qtrl optimizers, and keep the above linked functionality for goat/jopt? (which is why it was written i guess). and also verify how exactly qtrl is doing the time evolution, and why its a little bit different from paul's code

[ajgpitch]

I don't believe Qtrl ever had a concept of a density matrix as a specific target. Using the TraceDiff optimizer just works, but it does mean that the infidelity does match anything we might expect. It seems clear that the scale is wrong, but clearly this does not matter when simply looking to minimize it. I think this should be fixed in QOC in terms of what the Result reports. QOC can also check the inputs and set the correct scaling factor.

The evolution mismatch seems more serious and requires more investigation. I would suggest trying to find some minimum level example, perhaps with just one or two timeslots, to try and identify the cause.

[rochisha0]

Hi @nwlambert,

According to me the final_state functionality was originally implemented for the crab and grape optimizers. In contrast, the _goat optimizer already returns the final_state as part of its result object.

The PR referenced above correctly applies the scale_factor as intended. I cannot figure out how to return the grape final_state from qoc.

[ajgpitch]

The final state in Qtrl is saved by result.evo_full_final = final_evo on line 789 of optimizer.py. You can see there that this comes from a property of the dynamics object.
I think it would also make sense to set QOC's result.final_state to the same

[nwlambert]

just to complicate things in the QOC result object evo_full_final is a property, but it seems to just return final_states

qutip-qoc/src/qutip_qoc/result.py

Line 400 in d5ff028

def evo_full_final(self):

[rochisha0]

Yes, @nwlambert . I'm not sure why is qoc trying hard not to use direct results from qtrl.