g_a_integrals: possible missing a coefficient of lambda/S ?

[daviehh]

First, thanks for open sourcing this implementation! I was writing my own implementation in julia from scratch following your paper SciPostPhys.7.6.080, but noticed a discrepancy in the J values: basically, e.g. with this example the infidelity J_T agrees but the integral of g_a differs. Glancing over the code, look like the integral of g_a

krotov/src/krotov/optimize.py

Line 474 in 9f9a223

g_a_integrals[i_pulse] += abs(Δϵ) ** 2 * dt # dt may vary!

does not agree with the second line of Eq. 7 in the paper: the integrand missed the \lambda_{a,l}/S_l(t) term and is just (\delta \epsilon)^2, and when I also drop that \lambda_{a,l}/S_l(t) term my julia code does agree with your python code. If I understand it correctly, this does not affect the shape of the controls (may change a bit due to iterating process stopping at a different iteration) but may just be that the actual J values would be different from the definition in the reference paper? What would be the correct term for the running cost int g_a part of J?

Again, thanks a lot!

[goerz]

I was writing my own implementation in julia from scratch following your paper SciPostPhys.7.6.080

I got you covered there! Have a look at https://github.com/JuliaQuantumControl/Krotov.jl

It's not as polished yet as the Python version and still in "alpha" mode (no stable API), but it has some extra features such as parametrized pulses from the start, and of course huge gains in parallel performance.

Feel free to get in touch about the Julia implementation...

About the discrepancy, at first glance, I'd say you're right, but let me double-check. In any case, this only affects the printed output at each iteration, but does not change the optimization in any way. Let me check also what it means for ΔJ (what exactly Krotov guarantees to be smaller than zero for "monotonic convergence"). If you're right, we gotta be a bit careful about not dividing by zero (since S(t) is usually zero at the edges)

[daviehh]

@goerz Thanks for the quick reply! I was writing the julia code as part of a learining process. Just to add the λ/S term is also present in other refs, e.g. Eq. 4 in https://arxiv.org/pdf/1008.5126.pdf

I understand that this is should not change the pulse shapes and is more-or-less a printing output thing, but it may change the final shape a bit if the resulting J value is used as part of the logic in when to stop the iterations (haven't checked yet)?

Thank you for double-checking! Concerning the divide-by-zero if this is indeed an issue: it should indeed be dealt with carefully in the code, but not a mathematical issue since Δε is given by S/λ*Im[...], inserting it into g_a = λ/S*Δε^2 as follows would not give divide-by-zero when the code is written as the last line

[goerz]

Yes, that's all correct :-)

Do you want to have a look at #97 and compare with the numbers you're getting with your version?

[daviehh]

Thanks for the confirmation and fix! Now I get the same numbers for the J, J_T and ∫gₐ(t)dt with #97

One other thing I noticed is with the printout for the iterations is that the ΔJ column does not mean how much ΔJ changes between iterations for iteration after the first one, the code

krotov/src/krotov/info_hooks.py

Lines 586 to 591 in 9f9a223

	Σgₐdt = np.sum(kwargs['g_a_integrals'])
	J = J_T_val + Σgₐdt
	if iteration > 0:
	J_T_prev_val = J_T_prev(**kwargs)
	ΔJ_T = J_T_val - J_T_prev_val
	ΔJ = ΔJ_T + Σgₐdt

suggests that the ΔJ column is ΔJ_T plus ∫gₐ(t)dt and not difference of the J value between the iterations. Is this as intended (see below for example with 01_example_simple_state_to_state.ipynb)?

Thank you!

[goerz]

That's right, and very much intended: The point is that gₐ more generally includes a "reference field", which just "happens to be" the guess field for that iteration. However, to calculate a ΔJ, you have to use the same reference field in both the terms, so the Δϵ vanishes in one of them.

I could have sworn that was explained somewhere in the documentation, but I can't find it at the moment...

Anyway, thanks for checking! I'll merge this, then.

[daviehh]

Got it, thanks for explaining!