krotov.objectives module

Routines for formulating objectives.

Objectives, represented as an Objective instance, describe the physical objective of an optimization, e.g. a state-to-state transformation, or a quantum gate. This is distinct from the mathematical formulation of an optimization functional (krotov.functionals). For the same physical objective, there are usually several different functionals whose minimization achieve that objective.




A single objective for optimization with Krotov’s method.



Extend objectives for an “ensemble optimization”


Construct a list of objectives for optimizing towards a quantum gate


Convert Hamiltonian and Lindblad operators into a Liouvillian.

__all__: Objective, ensemble_objectives, gate_objectives, liouvillian


krotov.objectives.FIX_QUTIP_932 = False

Workaround for QuTiP issue 932.

If True, in Objective.mesolve(), replace any array controls with an equivalent function. This results in a signficant slowdown of the propagation, as it circumvents the use of Cython. Defaults to False on Linux, and True on any non-Linux system.

class krotov.objectives.Objective(*, initial_state, H, target, c_ops=None)[source]

Bases: object

A single objective for optimization with Krotov’s method.



>>> H0 = - 0.5 * qutip.operators.sigmaz()
>>> H1 = qutip.operators.sigmax()
>>> eps = lambda t, args: ampl0
>>> H = [H0, [H1, eps]]
>>> krotov.Objective(
...     initial_state=qutip.ket('0'), target=qutip.ket('1'), H=H
... )
Objective[|Ψ₀(2)⟩ to |Ψ₁(2)⟩ via [H₀[2,2], [H₁[2,2], u₁(t)]]]

ValueError – If any arguments have an invalid type or structure. This can be surpressed by setting the type_checking class attribute to False.


Giving collapse operators via c_ops only makes sense if the propagator passed to optimize_pulses() takes them into account explicitly. It is strongly recommended to set H as a Lindblad operator instead, see liouvillian().


The (time-dependent) Hamiltonian or Liouvillian in nested-list format, cf. qutip.mesolve.mesolve(). This includes the control fields.


qutip.Qobj or list


The initial state, as a Hilbert space state, or a density matrix.




An object describing the “target” of the optimization, for the dynamics starting from initial_state. Usually, this will be the target state (the state into which initial_state should evolve). More generally, it can be an arbitrary object meeting the conventions of a specific chi_constructor function that will be passed to optimize_pulses().


List of collapse operators, cf. mesolve(), in lieu of H being a Liouvillian.


list or None

str_use_unicode = True

Whether the string representation of an Objective may use unicode symbols, cf. summarize() (class attribute).

type_checking = True

By default, instantiating Objective with invalid types raises a ValueError. Setting this to False disables type checks in the initializer, allowing certain advanced use cases such as using plain numpy objects instead of QuTiP objects (class attribute).


The Objective containing the adjoint of all components.

This does not affect the controls in H: these are assumed to be real-valued. Also, will be left unchanged if its adjoint cannot be calculated (if it is not a target state).

mesolve(tlist, *, rho0=None, H=None, c_ops=None, e_ops=None, args=None, **kwargs)[source]

Run qutip.mesolve.mesolve() on the system of the objective.

Solve the dynamics for the H and c_ops of the objective, starting from the objective’s initial_state, by delegating to qutip.mesolve.mesolve(). Both the initial state and the dynamical generator for the propagation can be overridden by passing rho0 and H/c_ops.

  • tlist (numpy.ndarray) – array of time grid points on which the states are defined

  • rho0 (qutip.Qobj or None) – The initial state for the propagation. If None, the initial_state attribute is used.

  • H (qutip.Qobj or None) – The dynamical generator (Hamiltonian or Liouvillian) for the propagation. If None, the H attribute is used.

  • c_ops (list or None) – List of collapse (Lindblad) operators. If None, the c_ops attribute is used.

  • e_ops (list or None) – A list of operators whose expectation values to calculate, for every point in tlist. See qutip.mesolve.mesolve().

  • args (dict or None) – dictionary of parameters for time-dependent Hamiltonians and collapse operators

  • **kwargs – All further arguments will be passed to qutip.mesolve.mesolve().


Result of the propagation, see qutip.mesolve.mesolve() for details.

Return type


propagate(tlist, *, propagator, rho0=None, H=None, c_ops=None, e_ops=None, args=None, expect=<function expect>)[source]

Propagate the system of the objective over the entire time grid.

Solve the dynamics for the H and c_ops of the objective. If rho0 is not given, the initial_state will be propagated. This is similar to the mesolve() method, but instead of using qutip.mesolve.mesolve(), the propagate function is used to go between points on the time grid. This function is the same as what is passed to optimize_pulses(). The crucial difference between this and mesolve() is in the time discretization convention. While mesolve() uses piecewise-constant controls centered around the values in tlist (the control field switches in the middle between two points in tlist), propagate() uses piecewise-constant controls on the intervals of tlist (the control field switches on the points in tlist). The function expect is used to calculate expecation values; it receives two parameters, an operator from e_ops and a state, and must returnd the expectation value of the operator.

Comparing the result of mesolve() and propagate() allows to estimate the “time discretization error”. If the error is significant, a shorter time step shoud be used.


Result of the propagation, using the same structure as mesolve().

Return type


classmethod reset_symbol_counters()[source]

Reset the internal symbol counters used for printing objectives.

See summarize().

summarize(use_unicode=True, reset_symbol_counters=False)[source]

Return a one-line summary of the objective as a string.

  • use_unicode (bool) – If False, only use ascii symbols in the output

  • reset_symbol_counters (bool) – If True, reset the internal object counters (see reset_symbol_counters()) before calculating the result

The summarize() method (which is also used for the repr() and __str__ of an Objective) keeps per-process internal counters for the various categories of objects that may occur as attributes of an Objective (kets, bras, Hermitian operators, non-Hermitian Operators, density matrices, Liouvillians, Lindblad operators, numpy arrays, control functions). This allows to keep track of objects across multiple objectives. The counters can be reset with reset_symbol_counters().

The ouput uses various unicode symbols (or ascii-equivalents, if use_unicode is False):

  • ‘Ψ’ (‘Psi’) for qutip.Qobj quantum states (kets or bras)

  • ‘ρ’ (‘rho’) for qutip.Qobj operators that occur as initial or target states (density matrices)

  • ‘L’ for Lindblad operators (elements of c_ops)

  • ‘H’ for Hermitian qutip.Qobj operators (Hamiltonians)

  • ‘A’ for non-Hermitian qutip.Qobj operators in H

  • ‘𝓛’ (‘Lv’) for qutip.Qobj super-operators (Liouvillians)

  • ‘a’ for numpy arrays (of any dimension)

  • ‘u’ for (callable) control functions.


>>> from qutip import ket, tensor, sigmaz, sigmax, sigmap, identity
>>> u1 = lambda t, args: 1.0
>>> u2 = lambda t, args: 1.0
>>> a1 = np.random.random(100) + 1j*np.random.random(100)
>>> a2 = np.random.random(100) + 1j*np.random.random(100)
>>> H = [
...     tensor(sigmaz(), identity(2)) +
...     tensor(identity(2), sigmaz()),
...     [tensor(sigmax(), identity(2)), u1],
...     [tensor(identity(2), sigmax()), u2]]
>>> C1 = [[tensor(identity(2), sigmap()), a1]]
>>> C2 = [[tensor(sigmap(), identity(2)), a2]]
>>> ket00 = ket((0,0))
>>> ket11 = ket((1,1))
>>> obj = Objective(
...     initial_state=ket00,
...     target=ket11,
...     H=H
... )
>>> obj.reset_symbol_counters()
>>> obj.summarize()
'|Ψ₀(2⊗2)⟩ to |Ψ₁(2⊗2)⟩ via [H₀[2⊗2,2⊗2], [H₁[2⊗2,2⊗2], u₁(t)], [H₂[2⊗2,2⊗2], u₂(t)]]'
>>> obj = Objective(
...     initial_state=ket00,
...     target=ket11,
...     H=H,
...     c_ops=[C1, C2]
... )
>>> obj.summarize()
'|Ψ₀(2⊗2)⟩ to |Ψ₁(2⊗2)⟩ via {H:[H₀[2⊗2,2⊗2], [H₁[2⊗2,2⊗2], u₁(t)], [H₂[2⊗2,2⊗2], u₂(t)]], c_ops:([[L₀[2⊗2,2⊗2], a₀[100]]],[[L₁[2⊗2,2⊗2], a₁[100]]])}'
>>> obj.summarize(use_unicode=False)
'|Psi0(2*2)> to |Psi1(2*2)> via {H:[H0[2*2,2*2], [H1[2*2,2*2], u1(t)], [H2[2*2,2*2], u2(t)]], c_ops:([[L0[2*2,2*2], a0[100]]],[[L1[2*2,2*2], a1[100]]])}'
>>> copy.deepcopy(obj).summarize()  # different objects!
'|Ψ₂(2⊗2)⟩ to |Ψ₃(2⊗2)⟩ via {H:[H₃[2⊗2,2⊗2], [H₄[2⊗2,2⊗2], u₁(t)], [H₅[2⊗2,2⊗2], u₂(t)]], c_ops:([[L₂[2⊗2,2⊗2], a₂[100]]],[[L₃[2⊗2,2⊗2], a₃[100]]])}'
>>> copy.deepcopy(obj).summarize(reset_symbol_counters=True)
'|Ψ₀(2⊗2)⟩ to |Ψ₁(2⊗2)⟩ via {H:[H₀[2⊗2,2⊗2], [H₁[2⊗2,2⊗2], u₁(t)], [H₂[2⊗2,2⊗2], u₂(t)]], c_ops:([[L₀[2⊗2,2⊗2], a₀[100]]],[[L₁[2⊗2,2⊗2], a₁[100]]])}'
krotov.objectives.gate_objectives(basis_states, gate, H, *, c_ops=None, local_invariants=False, liouville_states_set=None, weights=None, normalize_weights=True)[source]

Construct a list of objectives for optimizing towards a quantum gate

  • basis_states (list[qutip.Qobj]) – A list of \(n\) canonical basis states

  • gate – The gate to optimize for, as a \(n \times n\) matrix-like object (must have a shape attribute, and be indexable by two indices). Alternatively, gate may be the string ‘perfect_entangler’ or ‘PE’, to indicate the optimization for an arbitrary two-qubit perfect entangler.

  • H (list or qutip.Qobj) – The Hamiltonian (or Liouvillian) for the time evolution, in nested-list format.

  • c_ops (list or None) – A list of collapse (Lindblad) operators, or None for unitary dynamics or if H is a Liouvillian (preferred!)

  • local_invariants (bool) – If True, initialize the objectives for an optimization towards a two-qubit gate that is “locally equivalent” to gate. That is, the result of the optimization should implement gate up to single-qubit operations.

  • liouville_states_set (None or str) – If not None, one of “full”, “3states”, “d+1”. This sets the objectives for a gate optimization in Liouville space, using the states defined in Goerz et al. New J. Phys. 16, 055012 (2014). See Examples for details.

  • weights (None or list) – If given as a list, weights for the different objectives. These will be added as a custom attribute to the respective Objective, and may be used by a particular functional (chi_constructor). The intended use case is for the liouville_states_set values ‘3states’, and ‘d+1’, where the different objectives have clear physical interpretations that might be given differing importance. A weight of 0 will completely drop the corresponding objective.

  • normalize_weights (bool) – If True, and if weights is given as a list of values, normalize the weights so that they sum to \(N\), the number of objectives. IF False, the weights will be used unchanged.


The objectives that define the optimization towards the gate. For a “normal” gate with a basis in Hilbert space, the objectives will have the basis_states as each initial_state and the result of applying gate to the basis_states as each target.

For an optimization towards a perfect-entangler, or for the local_invariants of the given gate, each initial_state will be the Bell basis state described in “Theorem 1” in Y. Makhlin, Quantum Inf. Process. 1, 243 (2002), derived from the canonical basis_states. The target will be the string ‘PE’ for a perfect-entanglers optimization, and gate for the local-invariants optimization.

Return type



ValueError – If gate, basis_states, and local_invariants are incompatible, or gate is invalid (not a recognized string)


The dimension of the basis_states is not required to be the dimension of the gate; the basis_states may define a logical subspace in a larger Hilbert space.


  • A single-qubit gate:

    >>> from qutip import ket, bra, tensor
    >>> from qutip import sigmaz, sigmax, sigmay, sigmam, identity
    >>> basis = [ket([0]), ket([1])]
    >>> gate = sigmay()  # = -i|0⟩⟨1| + i|1⟩⟨0|
    >>> H = [sigmaz(),[sigmax(), lambda t, args: 1.0]]
    >>> objectives = gate_objectives(basis, gate, H)
    >>> assert objectives == [
    ...     Objective(
    ...         initial_state=basis[0],
    ...         target=(1j * basis[1]),
    ...         H=H
    ...     ),
    ...     Objective(
    ...         initial_state=basis[1],
    ...         target=(-1j * basis[0]),
    ...         H=H
    ...     )
    ... ]
  • An arbitrary two-qubit perfect entangler:

    >>> basis = [ket(n) for n in [(0, 0), (0, 1), (1, 0), (1, 1)]]
    >>> H = [
    ...     tensor(sigmaz(), identity(2)) +
    ...     tensor(identity(2), sigmaz()),
    ...     [tensor(sigmax(), identity(2)), lambda t, args: 1.0],
    ...     [tensor(identity(2), sigmax()), lambda t, args: 1.0]]
    >>> objectives = gate_objectives(basis, 'PE', H)
    >>> from weylchamber import bell_basis
    >>> for i in range(4):
    ...     assert objectives[i] == Objective(
    ...        initial_state=bell_basis(basis)[i],
    ...        target='PE',
    ...        H=H
    ...     )
  • A two-qubit gate, up to single-qubit operation (“local invariants”):

    >>> objectives = gate_objectives(
    ...     basis, qutip.gates.cnot(), H, local_invariants=True
    ... )
    >>> for i in range(4):
    ...     assert objectives[i] == Objective(
    ...        initial_state=bell_basis(basis)[i],
    ...        target=qutip.gates.cnot(),
    ...        H=H
    ...     )
  • A two-qubit gate in a dissipative system tracked by 3 density matrices:

    >>> L = krotov.objectives.liouvillian(H, c_ops=[
    ...     tensor(sigmam(), identity(2)),
    ...     tensor(identity(2), sigmam())])
    >>> objectives = gate_objectives(
    ...     basis, qutip.gates.cnot(), L,
    ...     liouville_states_set='3states',
    ...     weights=[20, 1, 1]
    ... )

    The three states, for a system with a logical subspace of dimension \(d\) with a basis \(\{\ket{i}\}\), \(i \in [1, d]\) are:

    \[\begin{split}\Op{\rho}_1 &= \sum_{i=1}^{d} \frac{2 (d-i+1)}{d (d+1)} \ketbra{i}{i} \\ \Op{\rho}_2 &= \sum_{i,j=1}^{d} \frac{1}{d} \ketbra{i}{j} \\ \Op{\rho}_3 &= \sum_{i=1}^{d} \frac{1}{d} \ketbra{i}{i}\end{split}\]

    The explicit form of the three states in this example is:

    >>> assert np.allclose(objectives[0].initial_state.full(),
    ...     np.diag([0.4, 0.3, 0.2, 0.1]))
    >>> assert np.allclose(objectives[1].initial_state.full(),
    ...     np.full((4, 4), 1/4))
    >>> assert np.allclose(objectives[2].initial_state.full(),
    ...     np.diag([1/4, 1/4, 1/4, 1/4]))

    The objectives in this example are weighted (20/1/1):

    >>> "%.5f" % objectives[0].weight
    >>> "%.5f" % objectives[1].weight
    >>> "%.5f" % objectives[2].weight
    >>> sum_of_weights = sum([obj.weight for obj in objectives])
    >>> "%.1f" % sum_of_weights
  • A two-qubit gate in a dissipative system tracked by \(d + 1 = 5\) pure-state density matrices:

    >>> objectives = gate_objectives(
    ...     basis, qutip.gates.cnot(), L,
    ...     liouville_states_set='d+1'
    ... )

    The first four initial_states are the pure states corresponding to the Hilbert space basis

    >>> assert objectives[0].initial_state == qutip.ket2dm(ket('00'))
    >>> assert objectives[1].initial_state == qutip.ket2dm(ket('01'))
    >>> assert objectives[2].initial_state == qutip.ket2dm(ket('10'))
    >>> assert objectives[3].initial_state == qutip.ket2dm(ket('11'))

    The fifth state is \(\Op{\rho}_2\) from ‘3states’:

    >>> assert np.allclose(objectives[4].initial_state.full(),
    ...     np.full((4, 4), 1/4))
  • A two-qubit gate in a dissipative system tracked by the full Liouville space basis:

    >>> objectives = gate_objectives(
    ...     basis, qutip.gates.cnot(), L,
    ...     liouville_states_set='full'
    ... )

    The Liouville space basis states are all the possible dyadic products of the Hilbert space basis:

    >>> assert objectives[0].initial_state == ket('00') * bra('00')
    >>> assert objectives[1].initial_state == ket('00') * bra('01')
    >>> assert objectives[2].initial_state == ket('00') * bra('10')
    >>> assert objectives[3].initial_state == ket('00') * bra('11')
    >>> assert objectives[4].initial_state == ket('01') * bra('00')
    >>> assert objectives[5].initial_state == ket('01') * bra('01')
    >>> assert objectives[6].initial_state == ket('01') * bra('10')
    >>> assert objectives[7].initial_state == ket('01') * bra('11')
    >>> assert objectives[8].initial_state == ket('10') * bra('00')
    >>> assert objectives[9].initial_state == ket('10') * bra('01')
    >>> assert objectives[10].initial_state == ket('10') * bra('10')
    >>> assert objectives[11].initial_state == ket('10') * bra('11')
    >>> assert objectives[12].initial_state == ket('11') * bra('00')
    >>> assert objectives[13].initial_state == ket('11') * bra('01')
    >>> assert objectives[14].initial_state == ket('11') * bra('10')
    >>> assert objectives[15].initial_state == ket('11') * bra('11')
krotov.objectives.ensemble_objectives(objectives, Hs, *, keep_original_objectives=True)[source]

Extend objectives for an “ensemble optimization”

This creates a list of objectives for an optimization for robustness with respect to variations in some parameter of the Hamiltonian. The trick is to simply optimize over the average of multiple copies of the system (the Hs) sampling that variation. See Goerz, Halperin, Aytac, Koch, Whaley. Phys. Rev. A 90, 032329 (2014) for details.

  • objectives (list[Objective]) – The \(n\) original objectives

  • Hs (list) – List of \(m\) variations of the original Hamiltonian/Liouvillian

  • keep_original_objectives (bool) – If given as False, drop the original objectives from the result. This is especially useful if Hs contains the original Hamiltonian (which is often more straightforward)


List of \(n (m+1)\) new objectives that consists of the original objectives, plus one copy of the original objectives per element of Hs where the H attribute of each objectives is replaced by that element. Alternatively, for keep_original_objectives=False, list of \(n m\) new objectives without the original objectives.

Return type


krotov.objectives.liouvillian(H, c_ops)[source]

Convert Hamiltonian and Lindblad operators into a Liouvillian.

This is like qutip.superoperator.liouvillian(), but H may be a time-dependent Hamiltonian in nested-list format. H is assumed to contain a drift Hamiltonian, and the Lindblad operators in c_ops cannot be time-dependent.