qucontrol / krotov Goto Github PK

This issue is to open up the discussion about some details of specifying objectives:

• What to do when a target_state doesn't make sense (PE/LI optimization)
• How to handle initial states/target states not being in the same Hilbert space
• How to take into account transformations between lab and rotating frame when optimizing in the RWA
• Whether there are any control problems that cannot be formulated in the current implementation.

Currently, we specify the control objectives for Krotov's optimize_pulses by giving a list of Objective objects, where each Objective contains an initial_state, information about the dynamical generator (H, c_ops), and a target_state. This design mimics the QDYN implementation. Furthermore, it reflects that Krotov's method inherently can be parallelized across the propagation of different states (which we typically index by subscript k). The Python implementation already supports the parallelization of these independent "control trajectories".

For large classes of control problems, it is also natural to have a specific target_state associated with each specific "control trajectory". The assumption is that for every initial_state, the optimization goal is reached when the dynamical generator evolves it into the target_state. Examples include gate optimization, where the initial_states are the basis states and the target_states are the result of applying the gate to the basis states; and "robustness optimization", where multiple copies of the same initial_state/target_state are considered, under different dynamical generators.

As a side note, for gate optimization, there is a possible alternative that QuTiP's GRAPE implementation follows (probably due to its Matlab roots), where the "state" is the time evolution operator. That is, the initial_state is the identity and the target_state is the gate. Our Krotov implementation supports this setup as well, although the ability to parallelize is lost. In any case, it doesn't fundamentally alter anything.

There are some control problems where the use of a target_state is not natural. The only examples that I'm personally aware of are a gate optimization with a local-invariants functional, or the optimization towards a perfect entangler. Since we don't know in advance which locally-invariant gate or which perfect entangler the optimization will converge to, we cannot specify a target_state to which each basis state should be go.

I propose to allow the possibility to have the value None for the target_state attribute of an Objective. The target_states are not actually used for anything essential inside the Krotov implementation: It only uses the initial_states for the forward propagation. Then, the user-supplied chi_constructor routine that determines the boundary states for the backward propagation receives the entire list of objectives as an input, and can use it in any way it wants; the chi-routines for the standard gate-optimization functionals do look at the target_states defined in that list, but this is not a requirement. The chi_constructor routine can even bring along its own data separate from the objectives, e.g. via closures. This will be required for the local-invariants optimization, where the target gate will be stored inside the chi_constructor, not anywhere in the objectives.

The Python Krotov implementation currently uses the target_states for one thing: it calculates tau_vals as the complex overlap of the forward-propagated states and the target states. These values are passed to the info_hook and are also stored in the Result object. The motivation for calculating these is that for all optimizations where the definition of target_state makes sense, these quantities are extremely informative independent of the functional that is used, and indeed the standard functionals can usually directly be defined in terms of these tau_vals (thus saving work in the info_hook). However, the tau_vals are in no way essential to the algorithm, and we can easily handle the situation where the target_states are None by also setting the tau_vals to None.

In #4, @Basilewitsch made the point that it might be desirable to allow the initial_states and target_states to be defined in different Hilbert spaces, e.g. the initial_state on a bipartite system and the target_state in a subspace. This is currently prevented by the calculation of the tau_vals, which requires both states to be in the Hilbert space. My initial thought was that this is the way it should be: States are ultimate always in the full Hilbert space in which the dynamics take place, and since the chi_constructor can do whatever it wants, it is free to project/partially trace as appropriate. However, I also think that it is perfectly reasonable for a user to expect to be able to settarget_state in exactly the way @Basilewitsch was intending. I would propose to add an (optional) attribute transform to the Objective class to handle this. If set, the (callable) transform will be applied to the forward-propagated state before taking the overlap with the stored target_state. The transform could use QuTiP's eliminate_states method to change the shape of the state, which would solve the problem at hand.

The transform could also handle at least one other use case that we have in QDYN, but haven't considered in the Python implementation: the transformation between the rotating frame and the lab frame when optimizing in the rotating wave approximation. In this case, it is natural to want to write the initial_states and target_states in the lab frame, and to use the standard chi_constructor routines also written in the lab frame. For this to work, the forward-propagated states must be transformed from the rotating frame to the lab frame before calculating the chis.

So, while the transform is very strictly superfluous, it gives an elegant solution to at least these two problems, and possibly more, so I think we should include it.

Lastly, QDYN includes the possibility to use the functional J_op for "optimization of the expectation value of a sum of operators". I think this was added by @danielreich a long time ago. It seems to me that this is just another example of an optimization where target_states are not appropriate, but where the operator(s) whose expectation value should be maximized can be contained in the chi_constructor, just like the gate in the local-invariants optimization. Thus, there is no need to change the definition of Objective. Maybe @danielreich can confirm this?

In summary, my proposal at this point is the following:

• Allow target_state=None in Objective. The tau_vals in this case will also be None.
• Add an optional callable attribute transform to Objective that allows to apply an arbitrary transformation to the fw_states_T before passing them to chi_constructor. This allows for target_state to be in a different Hilbert space than initial_state, and it helps in working in the RWA.

Is is possible to use the "list string format" for time-dependent operators, and an input in optimize_pulses? If yes, note this in qutip_usage.rst and add an example.

In 3114718, I've rewritten the section "Other Optimization Methods" (formerly, "Krotov vs GRAPE and CRAB"). This includes separating out the comparison of Krotov vs GRAPE with respect to the method as such, and the somewhat limited implementation in QuTiP, as well as a more complete discussion of the numerical costs associated with gradient-free optimization.

@karl-horn, can you check whether your use of NLopt is covered sufficiently in the current version? Based on https://iopscience.iop.org/article/10.1088/1367-2630/aaf360/meta, it looked like you're using Subplex only as a direct alternative to Nelder-Mead. Please extend the documentation if there's something more fundamental to discuss about methods provided by the NLopt package, or if it should be included in the decision tree diagram in the last section in some way.

If a control is a callable that returns integers (most notably, 0), the pulse is internally allocated as an array of ints, instead of floats. The pulse updates are then also truncated to int values.

Consider these two notebooks illustrating the problem (thanks to @MatthiKrauss for reporting):

notebooks.zip

We need a Python implementation of QDYN's make_ensemble_targets: Take a list of objectives that all use the same Hamiltonian, and a list of Hamiltonians, and return a new list of objectives that duplicates the original objective once for each new Hamiltonian.

QuTiP version 4.4 was released recently and appears to break a number of krotov's tests.

See https://gist.github.com/goerz/ba8575b060497ec05d5a928925432b19 (run locally on a Macbook) or the Travis log at https://travis-ci.org/qucontrol/krotov/jobs/558664343.

The main problem seems to be that mesolve now prints extra output when it is used to propagate a Schrödinger equation (qutip/qutip#1047).

Add an example notebook that performs an ensemble optimization of a Rydberg gate, as in Phys. Rev. A 90, 032329 (2014)

This will be easier with #10

Like #3, but using density matrices, and spontaneous decay from level 2. This should yield solutions similar to analytical StiRAP.

At a technical level, this illustrates the use of density matrices.

This wasn't documented properly, but the intention was that the tau_vals that get passed to the chi_constructor are for the current iteration only. Instead, optimize_pulses currently is passing it the list of lists of tau_vals, for all iterations. In aa91f33 @FernandoGago and @MatthiKrauss adapted chis_ss and chis_sm to this wrong behavior. But the bug is really in optimize_pulses, and in the missing documentation what the interface for the chi_constructors is.

I will fix both of these.

Right now, custom tex macros like \ket in the documentation seem to work locally, but only after reloading the page multiple times, and they don't seem to work at all on RTD.

The perfect-entanglers needs very specific initial/target states

This is the equivalent to set_oct_LI_targets in QDYN.

See comment in qutip_usage.rst.

There should be a way to plug in optimized pulses into the objective, in order to propagate them easily.

This could also address qutip/qutip#932

Add an example notebook that optimizes the transfer from state 1 to state 3 in a Lambda-sytem, with a Pump and a Stokes laser, in the rotating wave approximation.

Dissipation from level 2 is modeled by a non-Hermitian Hamiltonian (decay term on level 2). Physically, this is the same as #4, but using a simpler numerical method. It should yield similar results as #4, and compare to the non-dissipative version #3.

At a technical level, this ensures that the backward propagation correctly uses the conjugate Hamiltonian.

The documentation needs to be synchronized with the current state of preprint (https://arxiv.org/abs/1902.11284v4)

The overlap routine in Krotov blindly assumes to be working on (Hermitian) density matrices. This was a deliberate choice (36b59aa), because calculating an unnecessary .dag() can be measurably expensive. However, for an optimization with the full Liouville basis the propagated states are not Hermitian, and a wrong overlap could cause serious problems. The Qobj.isherm flag could be used to decide whether to call dag() or not.

>>> import numpy as np
>>> import qutip
>>> from qutip import ket
>>> import krotov
>>> Qmagic = (1.0 / np.sqrt(2.0)) * qutip.Qobj(
...     np.array(
...         [[1,  0,  0,  1j], [0, 1j,  1,  0], [0, 1j, -1,  0], [1,  0,  0, -1j]],
...         dtype=np.complex128,
...     ),
...     dims=[[2, 2], [2, 2]],
... )
>>> def ketbra(a, b):
...     return ket(a) * ket(b).dag()
...
>>> rho_2 = ketbra('01', '10')
>>> (Qmagic.dag() * rho_2).tr()  # correct
0.7071067811865475
>>> (Qmagic * rho_2).tr()  # wrong
0.7071067811865475j
>>> krotov.second_order._overlap(Qmagic, rho_2)  # wrong
0.7071067811865475j

The check_convergence parameter that optimize_pulses should be able to receive to determine if the optimization can stop before iter_stop is reached is currently ignored.

In 262aea3, I've rewritten the section "Krotov's Method" in the documentation according to our discussion from earlier this week.

• @FernandoGago @karl-horn @MatthiKrauss, please review the section to check if you find it easier to understand now, especially in regard to the use of the second order, and the time discretization

• @danielreich Can you give this another quick review for mathematical correctness and consistent use of indices?

You may check off the above two boxes when you've finished the review, so that I can close this issue.

Please feel free to directly rewrite anything that you feel can be expressed more clearly, or any typos you see.

The current implementation of Krotov works for real pulses only. This needs to be fixed.

This is a prerequisite for #3, #4, #9 working as intended.

• Krotov version: dev version
• Python version: 3.6.8
• Operating System: macOS 10.13

Description

Currently the optimization method has no limit on the parameter space. For example, in the case of "simple state to state", the S(t) is initially set to be within [0, 1] but this range no longer holds. This agrees with practical consideration where the maximal B_x is limited.

Use the truncated transmon example from PRA 91, 062307 (2015), Eq (16).

This is the primary example for the second-order update.

Depends on #5, #6, #7

The quantity ∫gₐ(t)dt is a very useful indicator for judging the convergence, and whether λₐ was chosen appropriately.

The optimize_pulses routine should calculate this quantity and make it available to the info_hook. Also there should be an info_hook that mimics the output of the optimization in QDYN (show in example 1).

• Krotov version: 0.3.0
• Python version: 3.6.5
• Operating System: win10

Description

In mu.py
line 132,133

for i in ham_mapping[1:]:
    mu += (1j * eqm_factor) * objective.H[ham_mapping[i]][0]

ham_mapping[i] should be replaced by i

For convenience, it would be nice to have a mesolve method for the Objective class, to propagate the initial state, or any other given state.

There should also be a propagate routine that gets passed a propagator (just like optimize_pulses). This way, one could easily compare the two propagations to check for "discretization errors".

This is especially useful in combination with #16

Add a trivial example notebook for the state-to-state optimization in a two-level system

Ensure that all examples follow the recently corrected equation in How to optimize complex control fields

• 02_example_lambda_system_rwa_complex_pulse.ipynb @karl-horn
• 03_example_lambda_system_rwa_non_hermitian.ipynb @FernandoGago
• 06_example_3states.ipynb @goerz

Currently, the text representation of a list of objectives only number the controls. For example, for an ensemble optimization:

[Objective[|(3)⟩ - {[Herm[3,3], [Herm[3,3], u1(t)], [Herm[3,3], u2(t)], [Herm[3,3], u3(t)], [Herm[3,3], u4(t)]]} - |(3)⟩],
Objective[|(3)⟩ - {[Herm[3,3], [Herm[3,3], u1(t)], [Herm[3,3], u2(t)], [Herm[3,3], u3(t)], [Herm[3,3], u4(t)]]} - |(3)⟩],
Objective[|(3)⟩ - {[Herm[3,3], [Herm[3,3], u1(t)], [Herm[3,3], u2(t)], [Herm[3,3], u3(t)], [Herm[3,3], u4(t)]]} - |(3)⟩],
Objective[|(3)⟩ - {[Herm[3,3], [Herm[3,3], u1(t)], [Herm[3,3], u2(t)], [Herm[3,3], u3(t)], [Herm[3,3], u4(t)]]} - |(3)⟩],
Objective[|(3)⟩ - {[Herm[3,3], [Herm[3,3], u1(t)], [Herm[3,3], u2(t)], [Herm[3,3], u3(t)], [Herm[3,3], u4(t)]]} - |(3)⟩]]

It would be better if we also automatically numbered all the states and operators, e.g.

[Objective[|Ψ₀(3)⟩ to |Ψ₁(3)⟩ via [H₀[3,3], [H₁[3,3], u₀(t)], [H₂[3,3], u₁(t)], [H₃[3,3], u₂(t)], [H₄[3,3], u₃(t)]]],
 Objective[|Ψ₀(3)⟩ to |Ψ₁(3)⟩ via [H₀[3,3], [H₅[3,3], u₀(t)], [H₆[3,3], u₁(t)], [H₇[3,3], u₂(t)], [H₈[3,3], u₃(t)]]],
 Objective[|Ψ₀(3)⟩ to |Ψ₁(3)⟩ via [H₀[3,3], [H₉[3,3], u₀(t)], [H₁₀[3,3], u₁(t)], [H₁₁[3,3], u₂(t)], [H₁₂[3,3], u₃(t)]]],
 Objective[|Ψ₀(3)⟩ to |Ψ₁(3)⟩ via [H₀[3,3], [H₁₃[3,3], u₀(t)], [H₁₄[3,3], u₁(t)], [H₁₅[3,3], u₂(t)], [H₁₆[3,3], u₃(t)]]],
 Objective[|Ψ₀(3)⟩ to |Ψ₁(3)⟩ via [H₀[3,3], [H₁₇[3,3], u₀(t)], [H₁₈[3,3], u₁(t)], [H₁₉[3,3], u₂(t)], [H₂₀[3,3], u₃(t)]]]]

which would make it easy to verify that all objectives go from the same state Ψ₀ to the same state Ψ₁, and that the drift Hamiltonian H₀ is the same in all objectives.

There is also room for improvement in distinguishing density matrices from operators (ρ vs H (Herm)/A (NonHerm)), and making sure that target being an arbitrary object will never cause the printing to fail.

When calculating μ=∂H/∂ϵ for an ε that does not occur in H, the derivative_wrt_pulse function currently returns the integer 0 for μ, which is not compatible with how μ is used in optimize_pulses: we need μ to be a callable. In the zero-case, μ(state) must return a zero-state while preserving the dimensions of state.

In addition to value_below, there should also be a value_above in krotov.convergence. The documentation should make it clear that evaluating functionals (errors or fidelities) and deciding convergence based on the calculated values is not inherent to the optimization, but entirely up to whatever convention the user chooses via the info_hook the pass.

Also: de-emphasize the use of glom (too technical).

By using the Hilbert-Schmidt distance as state overlap measure for a gate optimization employing the 3-states method, the optimization is not monotonically minimizing the functional but keeps fluctuating. Interestingly, if the co-states chis[k] are exchanged by -chis[k], the optimization indeed shows monotonic convergence in terms of a maximization of the Hilbert-Schmidt distance. Moreover, the optimization minimizes the functional correctly, if the real-part functional is used.

The implementations of the Hilbert-Schmidt distance and its co-states construction seem correct to me. Maybe a problem with the (non-?)convexity of the Hilbert-Schmidt distance?

The current implementation of Krotov only allows for propagation of Hilbert space kets in the optimization routine. In order to proceed to dissipative systems, propagation with density matrices and the Liouville-von Neumann equation needs to be added. Note that 'mesolve' routine from QuTiP already allows that, i.e., propagation outside the Krotov optimization is possible.

This is a prerequisite for #4, #12.

Specially when talking about bounded variables, using an optimization for a Hamiltonian H( f ( Ɛ(t) )) can be really useful, with Ɛ(t) being the parameter to optimize and f a function with a bounded image.

So as Christiane and Daniel R. used in one of their past works, for a variable α with bounds [a,b] one could write α = b*( tanh(Ɛ) + 1)/2 + a and try to optimize Ɛ. Due to the monotonic behaviour in tanh I don't think second order Krotov would be necessary for this to work, so it would be a nice feature to have even in first order Krotov.

We need a complete set of example notebooks that covers all use cases of the package.
See https://krotov.readthedocs.io/en/latest/contributing.html#contribute-examples for how to add notebooks.

Based on today's discussion, we want at least the following examples:

• Perfect entangler optimization #8
• Optimization with complex pulses #3
• State-to-state density matrix optimization #4
• Ensemble optimization for robustness #11
• "3-states" optimization of dissipative quantum gate #12

We will skip the example for a state-dependent constraint, as this functionality will not be added for version 1.0.

This is a meta-issue tracking the open issues for each of the example notebooks. See also the associated project: https://github.com/qucontrol/krotov/projects/1

Add an example notebook for the optimization of a dissipative quantum gate using the "3 states" of New J. Phys. 16, 055012 (2014)

There are a few technical issues with Objectives:

• It appears copy.deepcopy taps into the pickling mechanism when __deepcopy__ is not explicitly defined, which results in callable controls being overwritten by placeholders during the copy.
• custom attributes (such as weight) are lost when copying an objective.
• obj1 == obj2 throws a ValueError if any of the controls are numpy arrays (since numpy arrays cannot be compared with ==).

For publication, we should go through the entire documentation carefully, checking for completeness, errors, and readability.

Independent of whether or not we choose to submit to JOSS, their checklist regarding a statement of purpose sounds useful:

• A list of the authors of the software
• Author affiliations
• A summary describing the high-level functionality and purpose of the software for a diverse, non-specialist audience
• A clear statement of need that illustrates the purpose of the software
• Mentions (if applicable) of any ongoing research projects using the software or recent scholarly publications enabled by it

Beyond that, I am planning a few more sections in the Howto's:

• How to optimize in a dissipative system (use Liouvillian, not c_ops)
• How to optimize for robust pulses (ensemble optimization)
• How to parallelize the optimization, cf #29
• How to maximize numerical efficiency (use Cython, show in #12)
• How to deal with the optimization running out of memory (out-of-memory arrays like dask)

The "X-Gate for a Transmon Qubit" example is a very simple example of a gate optimization, and thus it fits in quite well into into the progression of example notebooks, as the first example for a gate optimization. It is preceded by state-to-state optimizations and followed up by more complex gate optimization (dissipative "3 states" optimization, and perfect entanglers).

However, the example could need a little bit of cleanup, and an extended discussion.

Prerequisite for #13

Any math in the markdown of notebook gets displayed incorrectly in the Sphinx documentation, e.g. as

This is from the example notebook for #19 (commit f2ec03e)

• Krotov version: stable
• Python version: 3.7
• Operating System: MacOS High sierra

Description

I wanted to download a PDF version of the Documentation. But I couldn't find an option on the Readthe Docs site. I think you need to build the documentation at least once for me to download it. A great format for issues by the way 😃.

What I Did

Tried to download PDF of Documentation. What I expected
option on Read thedocs own Documentation site
What I found,
Option on krotov Documentation site

Sorry for the silly issue. I just have a problem with staring at the computer for too long. This would happen when reading the documentation on the computer.

Taking the adjoint of an objective that contains a numpy array as a control throws an AttributeError ('numpy.ndarray' object has no attribute 'dag')

There is an incongruity between the current definition of τ = ⟨ϕ(T)|ϕtgt⟩ and the definition of e.g. χss = τ ϕtgt and χsm = (∑τ) ϕtgt. As a result, Example 1 (stage-to-state transition in a TLS) does not converge monotonically if J_T_ss is used instead of J_T_re. It would work with χss = τ⁺ ϕtgt, but the better solution is to redefine τ as ⟨ϕtgt | ϕ(T)⟩. This is in fact the original definition from PalaoPRA2003.

Incidentally, QDYN's documentation is also confused about this: the τ are correctly calculated as ⟨ϕtgt | ϕ(T)⟩, but the documentation claims that they are ⟨ϕ(T)|ϕtgt⟩, in oct_funcs.F90. @Basilewitsch, can you make sure this gets added to QDYN's TODO list?

For the sake of completeness, we need F_sm, F_ss, J_T_sm, J_T_ss, chis_sm, chis_ss.

@MatthiKrauss, thanks for adding these in b8dea5c

Could you also add basic tests for these functionals in tests/test_functionals.py (similar to the existing tests for the other functionals)? The added routines should have complete coverage, see make coverage locally, or https://coveralls.io/github/qucontrol/krotov?branch=master after pushing.

Lastly, can you add the names of the new functions to the __all__ list in src/krotov/functionals.py? This is just the convention of declaring the "public" routines in a module.

• Implement the second order update equation
• Add code for numerical estimation of σ(t)

This should mainly point out that you can do spectral constraints simply by applying a spectral filter after each iteration, but if @danielreich want's to add a sentence about the "proper" spectral constraints, that could go in there as well (with the the caveat that this is a fairly big change in Krotov's method, and not supported or planned by the implementation in the package)

Krotov version: stable
Python version: 3.7
Operating System: Ubuntu 18.04

Description

I was trying to convert the example Jupyter notebooks
Describe what you were trying to get done.
Tell us what happened, what went wrong, and what you expected to happen.

What I Did

Paste the command(s) you ran and the output.
If there was a crash, please include the traceback here.

The CONTRIBUTING docs should explain how to write commit messages that automatically get associated with issues on Github. Doing this consistently is a great way to stay organized.

Currently, the guess pulses are calculated in such a way that the values at 0 and T are not preserved. This runs counter to the point in the documentation (http://nbviewer.jupyter.org/github/qucontrol/krotov/blob/4dff23f/docs/time_discretization.ipynb, input cell 6) that the guess controls must be discretized to tlist before converting them to guess pulses.

Calling optimize_pulses with parallel_map=qutip.parallel.parallel_map, e.g. in the example for the gate optimization (#27) does not result in the expected speedup by a factor of two. Instead, it appears to run much much slower than the non-parallelized version. This could be a bug in qutip, but we'll have to investigate. In any case, we need to make sure there is a way to run the optimization with parallelization, and the gate optimization seems like a good example where to show this (as well as the ensemble optimization example, #11)

The optimize_pulses function should have a continue_from option that can receive a Result object from a previous call to optimize_pulses, and continue the optimization from that previous result (presumably, with a larger value of iter_stop).

Furthermore, there should be an option make_snapshots=<n: int> that dumps the current Result object to disk every n iterations. This is to make Krotov robust with respect to crashes (or being killed by a shutdown/HPC scheduler): we'll be able to continue from the last dumped Result.

Add an example notebook that optimizes the transfer from state 1 to state 3 in a Lambda-sytem, with a Pump and a Stokes laser, without dissipation, and in the rotating wave approximation.

At a technical level, this is mainly to illustrate (and test) the use of complex control fields.

The current info_hook does more than just calculate values that end up in Result.info_vals: It is also the proper place to dynamically modify the objectives or the λₐ. While an advanced use case, this is very much intended, and the only way to to implement things like the following:

• Adjusting λₐ dynamically
• A "moving grid" Hamiltonian
• Filters that act on more than one pulse at the same time (e.g., when there is one pulse for the real part and one pulse for the imaginary part of a complex control, but a filter needs to be applied at the level of the complex control)

Generally speaking, the info_hook is where a user can inject arbitrary code into the optimization, after each iteration. To emphasize that the info_hook is not just for "information", it should be renamed to iter_hook.