.. _fermionic_synthesis_explanation:

Synthesize an Evolution gate in fermionic space
===============================================

.. important::

   The concepts in this guide are currently available only in the Python API.
   Equivalent functionality will be made available in the C API in a future release.

An :class:`.Evolution` gate carries the time evolution :math:`e^{-itH}` of a whole operator. A
`fermion-to-fermion` synthesis rewrites that gate into a :class:`.FermionicCircuit` of smaller
fermionic gates. Both sides of the rewrite stay in fermionic space, which is what distinguishes it
from the fermion-to-qubit synthesis that follows: the operator's
:attr:`~qiskit_fermions.operators.OperatorTrait.groups` (see :ref:`grouping_explanation`) survive
it, whereas a fermion-to-qubit mapping returns a
:class:`~qiskit.quantum_info.SparseObservable`, which has no concept of grouping.

Two knobs select the method. :attr:`.Evolution.synthesis` sets it on one gate, and the
:class:`.FermionicTrotterization` transpiler pass sets it across a whole circuit. Both take a
:class:`.FermionicEvolutionSynthesis`, the interface any such rewrite implements. The methods this
package currently provides are product formulas, which is what the rest of this guide uses:
:class:`.FermionicLieTrotter` for a first-order formula (the default) and
:class:`.FermionicSuzukiTrotter` for a higher even order. Anything else implementing the interface
is selected the same way, and the sections on the two knobs apply to it unchanged.

.. note::
   This fermion-to-fermion step is **optional**. An :class:`.Evolution` gate can go straight to the
   fermion-to-qubit stage, which maps and synthesizes it whole however many terms it holds.
   Splitting it first is a choice, made because the resulting gates are individually cheaper or
   because the split exposes structure that a later stage can use.

   The two steps compose, so an evolution can end up approximated at either level, at both, or at
   neither. A product formula approximates unless its factors commute; a synthesis method need not
   approximate at all.

.. note::
   Skipping this step does not have to mean losing the grouping entirely. Wrapping the mapper in
   :func:`.group_wise` maps the operator one group at a time and sums the results, so the mapped
   Pauli terms arrive grouped rather than interleaved. A qubit-side product formula synthesizes them
   in the order it receives, so commuting terms that end up adjacent can be scheduled together.

   On the Hamiltonian below that recovers the whole two-qubit depth benefit of splitting: 18 instead
   of 30, at an unchanged gate count. The difference is *where* the structure is expressed.
   Splitting in fermionic space makes each group a separate gate, which every later stage can see
   and act on; ``group_wise`` leaves one gate and only reorders the terms inside it, which is enough
   for scheduling but leaves nothing for a fermionic pass to work with. Use it when the evolution
   has to stay a single gate, and reach for a fermionic product formula otherwise.

Build a grouped Hamiltonian
---------------------------

Both product formulas split an operator into factors, and for a grouped operator those factors are
its groups. The grouping is therefore what determines both the cost and the accuracy of the result.
Take a six-site Fermi-Hubbard chain, grouped the way its terms commute: the hopping terms on even
bonds, the hopping terms on odd bonds, and the diagonal interaction.

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> from qiskit_fermions.operators import FermionOperator, cre, ann
   >>>
   >>> def fermi_hubbard_chain(num_sites, tunneling, interaction):
   ...     """A 1D Fermi-Hubbard chain grouped into even bonds, odd bonds, and the interaction."""
   ...     terms, groups = [], []
   ...     for offset in (0, 1):
   ...         for j in range(offset, num_sites - 1, 2):
   ...             terms += [
   ...                 ([cre(j), ann(j + 1)], -tunneling),
   ...                 ([cre(j + 1), ann(j)], -tunneling),
   ...             ]
   ...             groups += [offset, offset]
   ...     for j in range(num_sites - 1):
   ...         terms.append(([cre(j), ann(j), cre(j + 1), ann(j + 1)], interaction))
   ...         groups.append(2)
   ...     operator = FermionOperator.from_terms(terms)
   ...     operator.groups = groups
   ...     return operator
   >>>
   >>> num_sites = 6
   >>> hamiltonian = fermi_hubbard_chain(num_sites, tunneling=1.0, interaction=2.0)
   >>> hamiltonian.num_groups()
   3

Fifteen terms in three groups. Each hopping group pairs every term with its conjugate partner, which
matters: a factor has to be Hermitian for its exponential to be unitary, and that does not follow
from the sum being Hermitian. :func:`.groups_are_hermitian` checks that convention rather than
leaving it to inspection:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> from qiskit_fermions.operators.terms.grouping import groups_are_hermitian
   >>>
   >>> groups_are_hermitian(hamiltonian)
   [True, True, True]

Separating :math:`a^\dagger_0 a_1` from its partner :math:`a^\dagger_1 a_0` would break it, and
leaving the operator ungrouped has the same effect, since both formulas then split it *term by term*
and a lone :math:`a^\dagger_i a_j` is not Hermitian:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> ungrouped = fermi_hubbard_chain(num_sites, tunneling=1.0, interaction=2.0)
   >>> ungrouped.groups = range(len(list(ungrouped.iter_terms())))  # one group per term
   >>> all(groups_are_hermitian(ungrouped))
   False

Such a factor is a valid fermionic operator, and this package will build the gate for it, so nothing
complains until the mapped operator reaches
:class:`~qiskit.circuit.library.PauliEvolutionGate`, which requires real coefficients and raises
``ValueError: Operator contains complex coefficients, which are not supported``. Checking the
grouping up front turns that late failure into an answer available before any circuit is built.

Set the method on one gate
--------------------------

Pass a synthesis method to the gate, then decompose it. Each gate it emits comes back as its own
:class:`.Evolution` gate:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> from qiskit_fermions.circuit import FermionicCircuit
   >>> from qiskit_fermions.circuit.library import Evolution
   >>> from qiskit_fermions.circuit.library.synthesis import (
   ...     FermionicLieTrotter,
   ...     FermionicSuzukiTrotter,
   ... )
   >>>
   >>> total_time = 1.0
   >>>
   >>> def factors(synthesis):
   ...     circuit = FermionicCircuit(num_sites)
   ...     circuit.append(
   ...         Evolution(num_sites, hamiltonian, time=total_time, synthesis=synthesis),
   ...         circuit.modes,
   ...     )
   ...     return circuit.decompose().count_ops()["Evolution"]
   >>>
   >>> factors(FermionicLieTrotter())
   3
   >>> factors(FermionicSuzukiTrotter(order=2))
   5
   >>> factors(FermionicSuzukiTrotter(order=2, reps=4))
   20

First order emits one factor per group. Order two emits the symmetric palindrome, so the two outer
groups appear twice at half the time each, giving five. Four repetitions of that give twenty. Every
order splits the operator the same way; what changes is the ordering of the factors and their
evolution times, not the partition.

.. note::
   The factors are marked :attr:`~.Evolution.atomic`, so decomposing again leaves them in place
   rather than splitting them further. The decomposition has a fixed point, which is what makes it
   safe to call :meth:`~.FermionicCircuit.decompose` repeatedly.

Set the method for a whole pipeline
-----------------------------------

Setting :attr:`~.Evolution.synthesis` per gate means threading the choice through everything that
builds one, including :class:`.UCC` and :class:`.UCJ`, which construct their gates internally. The
:class:`.FermionicTrotterization` pass applies one method to every :class:`.Evolution` in a circuit
instead, so the choice is made once for the pipeline:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> from qiskit_fermions.transpiler import FermionicPassManager
   >>> from qiskit_fermions.transpiler.passes import FermionicTrotterization
   >>>
   >>> circuit = FermionicCircuit(num_sites)
   >>> circuit.append(Evolution(num_sites, hamiltonian, time=total_time), circuit.modes)
   >>>
   >>> expanded = FermionicPassManager(
   ...     FermionicTrotterization(FermionicSuzukiTrotter(order=2, reps=4))
   ... ).run(circuit)
   >>> expanded.count_ops()["Evolution"]
   20

Twenty gates: five factors per order-two sweep, times four repetitions -- the same count
:meth:`~.FermionicCircuit.decompose` gives for that method above. The pass carries the formula out
itself, so no separate expansion step is needed, and running it twice changes nothing, since the
factors it emits are :attr:`~.Evolution.atomic` and so carry no definition of their own.

.. note::
   Pass ``apply=False`` to only *select* the method and leave the expansion to something else, such
   as Qiskit's :class:`~qiskit.transpiler.passes.Decompose` (naming the gate as
   ``Decompose("Evolution")`` restricts it to the evolutions, leaving the other fermionic gates to
   their own synthesis plugins).

   Be aware of what an unexpanded gate does, though: it reaches the fermion-to-qubit stage whole,
   where it is mapped without :attr:`~.Evolution.synthesis` being consulted. That makes for an easy
   mistake when comparing methods, because every one of them then produces identical output and the
   comparison silently measures nothing. On the Hamiltonian above, leaving the order-two gate
   unexpanded turns the 200-gate circuit below into the same 30-gate circuit that first order gives.

Measure the payoff
------------------

A higher order costs depth, so the question is not whether it is more accurate but whether it is
more accurate *per gate*. Two numbers are needed per method: the error of the state it produces, and
the two-qubit cost of the circuit that produced it.

The error is measured with `ffsim`_, which simulates a :class:`.FermionicCircuit` **directly** (see
:ref:`ffsim_relationship_explanation`). No fermion-to-qubit mapping is involved, so the comparison
isolates the fermionic product formula from anything the mapping or the qubit-side synthesis might
contribute. It also works in the fixed-particle-number sector, which is both smaller than the full
:math:`2^N` space and the space the evolution actually stays in.

.. invisible-code-block: python

   >>> from qiskit_fermions.utils.optionals import HAS_FFSIM

.. skip: start if(not HAS_FFSIM)

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> import ffsim
   >>> import numpy as np
   >>> from scipy.sparse.linalg import expm_multiply
   >>> from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
   >>>
   >>> norb, nelec = num_sites, (3, 0)  # three spinless fermions on six modes
   >>> initial = ffsim.hartree_fock_state(norb, nelec)
   >>>
   >>> hamiltonian_operator = ffsim.linear_operator(hamiltonian, norb=norb, nelec=nelec)
   >>> exact = expm_multiply(
   ...     -1j * total_time * hamiltonian_operator,
   ...     initial,
   ...     traceA=-1j * total_time * ffsim.trace(hamiltonian, norb=norb, nelec=nelec),
   ... )
   >>>
   >>> def cost_and_error(synthesis):
   ...     circuit = FermionicCircuit(num_sites)
   ...     circuit.append(
   ...         Evolution(num_sites, hamiltonian, time=total_time, synthesis=synthesis),
   ...         circuit.modes,
   ...     )
   ...     # the method is set on the gate, so decompose() is what expands it here; a pipeline
   ...     # would reach for FermionicTrotterization instead (see the note above)
   ...     circuit = circuit.decompose()
   ...     evolved = ffsim.apply_unitary(initial, circuit, norb=norb, nelec=nelec)
   ...     error = 1 - abs(np.vdot(evolved, exact)) ** 2
   ...     qubit_circuit = generate_preset_jw_pass_manager().run(circuit).decompose(reps=6)
   ...     two_qubit = lambda instruction: len(instruction.qubits) == 2
   ...     return error, qubit_circuit.count_ops()["cx"], qubit_circuit.depth(two_qubit)
   >>>
   >>> methods = [
   ...     ("order 1, reps 1 ", FermionicLieTrotter()),
   ...     ("order 2, reps 1 ", FermionicSuzukiTrotter(order=2)),
   ...     ("order 1, reps 10", FermionicLieTrotter(reps=10)),
   ...     ("order 2, reps 4 ", FermionicSuzukiTrotter(order=2, reps=4)),
   ... ]
   >>> for label, synthesis in methods:
   ...     error, gates, depth = cost_and_error(synthesis)
   ...     print(f"{label}  error {error:.6f}  ({gates:3d} CX, depth {depth:3d})")
   order 1, reps 1   error 0.515997  ( 30 CX, depth  18)
   order 2, reps 1   error 0.180958  ( 50 CX, depth  26)
   order 1, reps 10  error 0.004026  (300 CX, depth 162)
   order 2, reps 4   error 0.000421  (200 CX, depth 104)

.. skip: end

Read the last two rows together, because they are the comparison that matters. Order two with four
repetitions is roughly ten times more accurate than order one with ten, and it gets there with a
third fewer entangling gates. Spending a gate budget on a higher order beats spending it on more
first-order steps here.

That result is not universal. Increasing :attr:`~.FermionicSuzukiTrotter.reps` and increasing
:attr:`~.FermionicSuzukiTrotter.order` reduce the error at different rates for the same kind of
cost, and which wins depends on the operator and on the time being evolved. The two rows above are
how to find out for a given problem.

.. note::
   The fermionic formula composes with the one chosen for the fermion-to-qubit stage (see
   :attr:`.MapperFnEvolutionSynthesis.product_formula`). Where both approximate, the accuracy of the
   result is governed by the weaker of the two, so raising the order here while the qubit-side
   formula stays first-order buys little.

Restrict the pass to the gates that benefit
-------------------------------------------

A higher order only helps when the factors do not commute. If they all commute, the first-order
formula is already exact and the extra factors are pure overhead. A purely diagonal operator makes
that concrete, since number operators commute with each other:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
   >>>
   >>> diagonal = FermionOperator.from_terms(
   ...     [([cre(j), ann(j), cre(j + 1), ann(j + 1)], 1.0 + 0.1 * j) for j in range(num_sites - 1)]
   ... )
   >>> diagonal.groups = list(range(num_sites - 1))
   >>>
   >>> def diagonal_cost(synthesis):
   ...     circuit = FermionicCircuit(num_sites)
   ...     circuit.append(
   ...         Evolution(num_sites, diagonal, time=total_time, synthesis=synthesis), circuit.modes
   ...     )
   ...     transpiled = generate_preset_jw_pass_manager().run(circuit.decompose())
   ...     return transpiled.decompose(reps=6).count_ops()["cx"]
   >>>
   >>> diagonal_cost(FermionicLieTrotter())
   10
   >>> diagonal_cost(FermionicSuzukiTrotter(order=2))
   18

Both circuits are exact, so the second spends 80% more entangling gates for nothing. The diagonal
Coulomb operators of a :class:`.UCJ` have this shape, which is why
:attr:`.FermionicTrotterization.filter` exists: it takes a predicate over the
:class:`~qiskit.dagcircuit.DAGOpNode` and leaves the rejected gates untouched, so one pass can raise
the order on the gates that benefit and leave the rest at first order.

A circuit holding both kinds of evolution shows what that saves. The predicate below accepts a gate
only when its operator moves particles between modes, which is exactly when the groups fail to
commute and a higher order has something to buy:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> def is_diagonal(operator):
   ...     """Whether every term only counts particles, rather than moving them between modes."""
   ...     return all(
   ...         sorted(mode for created, mode in actions if created)
   ...         == sorted(mode for created, mode in actions if not created)
   ...         for actions, _ in operator.iter_terms()
   ...     )
   >>>
   >>> def mixed_circuit():
   ...     circuit = FermionicCircuit(num_sites)
   ...     circuit.append(Evolution(num_sites, hamiltonian, time=total_time), circuit.modes)
   ...     circuit.append(Evolution(num_sites, diagonal, time=total_time), circuit.modes)
   ...     return circuit
   >>>
   >>> def cx_count(circuit):
   ...     return generate_preset_jw_pass_manager().run(circuit).decompose(reps=6).count_ops()["cx"]
   >>>
   >>> order_2 = FermionicSuzukiTrotter(order=2)
   >>>
   >>> cx_count(mixed_circuit())  # first order everywhere
   40
   >>> cx_count(FermionicPassManager(FermionicTrotterization(order_2)).run(mixed_circuit()))
   68
   >>> cx_count(
   ...     FermionicPassManager(
   ...         FermionicTrotterization(
   ...             order_2, filter=lambda node: not is_diagonal(node.op.operator)
   ...         )
   ...     ).run(mixed_circuit())
   ... )
   60

Raising the order on both gates costs 28 extra entangling gates; restricting it to the Hamiltonian
costs 20 and buys the same accuracy, because the eight the filter saved were spent on an evolution
that was already exact at first order.

.. caution::
   Neither knob verifies that the factors are Hermitian, and an incorrectly grouped operator can
   look *better* at an even order than at first order while still being wrong, because the
   symmetrization partially cancels the error. Group conjugate partners together, as
   ``fermi_hubbard_chain`` above does.

Next steps
----------

- Read :ref:`grouping_explanation` for how the groups these product formulas split on are assigned.
- Read :ref:`transpilation_explanation` for the stages a :class:`.FermionicCircuit` passes through.
- See :ref:`1d_fermi_hubbard` for a worked case where the fermionic formula supplies the step count
  that a custom qubit-side synthesis cannot.
- Consult :mod:`~qiskit_fermions.circuit.library.synthesis` for the synthesis methods themselves,
  and :class:`.QDriftTrotterization` for a randomized alternative.

.. _ffsim: https://qiskit-community.github.io/ffsim/