Sintetizar una puerta de evolución en el espacio bosónico

La evolución temporal $e^{-itH}$ de un hamiltoniano que es suma de piezas que no conmutan no se puede implementar exactamente puerta por puerta. Una fórmula de producto la aproxima mediante una secuencia de evoluciones cortas bajo cada una de las piezas, y el error de la aproximación se denomina error de Trotter.

Puedes aplicar una fórmula de producto en dos lugares, y no son intercambiables:

1. En el espacio bosónico, antes de haber elegido ningún truncamiento ni codificación. La fórmula ve los grupos del operador (consulta Agrupar términos del operador: usar la estructura del operador) y puede entrelazarlos.

2. En el espacio de qubits, dentro de PauliEvolutionGate de Qiskit, después de mapear el operador. La fórmula ve un operador de Pauli por grupo y solo puede subdividir los términos dentro de él.

Solo la primera puede reducir el error entre grupos que no conmutan, que suele ser el dominante. Esta guía muestra cómo usarla.

Importante

Ambos pasos son opcionales. Una puerta Evolution que nunca se descompone en el espacio bosónico llega entera a la etapa de bosón a qubit y se mapea y sintetiza allí. Eso no introduce error a nivel bosónico, solo el que introduzca la fórmula a nivel de qubits.

Preparación

Una cadena de Bose-Hubbard de dos sitios, el ejemplo más pequeño en el que los grupos no conmutan (consulta Simular la dinámica de Bose-Hubbard con una codificación a elegir para el modelo completo). Se agrupa en un grupo hermítico para el enlace y un grupo diagonal para la interacción en el sitio:

>>> import numpy as np
>>> from qiskit.quantum_info import SparsePauliOp, Statevector
>>> from qiskit.transpiler import PassManager
>>> from scipy.sparse import csc_matrix
>>> from scipy.sparse.linalg import expm_multiply
>>> from qiskit_bosons.circuit import BosonicCircuit
>>> from qiskit_bosons.circuit.library import Evolution, InitializeModes
>>> from qiskit_bosons.mappers import LevelSpec
>>> from qiskit_bosons.mappers.library import occupation_encoding
>>> from qiskit_bosons.operators import BosonOperator, ann, cre
>>> from qiskit_bosons.transpiler.presets import generate_preset_occupation_pass_manager
>>>
>>> num_sites, max_occupation, total_time = 2, 2, 0.8
>>> bonds = [(0, 1)]
>>> terms = {}
>>> for left, right in bonds:
...     terms[(cre(left), ann(right))] = -1.0
...     terms[(cre(right), ann(left))] = -1.0
>>> for site in range(num_sites):
...     terms[(cre(site), cre(site), ann(site), ann(site))] = 1.0
>>> hamiltonian = BosonOperator.from_dict(terms)
>>> hamiltonian.groups = [
...     len(bonds)
...     if len({mode for _, mode in term}) == 1
...     else bonds.index(tuple(sorted({mode for _, mode in term})))
...     for term, _ in hamiltonian.iter_terms()
... ]
>>> hamiltonian.num_groups()
2

Para medir un error necesitas una referencia. La referencia adecuada es la evolución exacta del hamiltoniano tal como se mapea con el mismo truncamiento, porque el truncamiento es una aproximación aparte:

>>> spec = LevelSpec(max_occupation, "unary_projected")
>>> mapped = occupation_encoding(hamiltonian, spec).simplify()
>>> matrix = csc_matrix(SparsePauliOp.from_sparse_observable(mapped).to_matrix(sparse=True))
>>> pm = generate_preset_occupation_pass_manager(spec)
>>> pm.qubit = PassManager()  # keep each group as one exact PauliEvolutionGate
>>>
>>> def build(**evolution_kwargs):
...     circuit = BosonicCircuit(num_sites)
...     circuit.append(InitializeModes([1, 1]), circuit.modes)
...     circuit.append(
...         Evolution(num_sites, hamiltonian, time=total_time, **evolution_kwargs), circuit.modes
...     )
...     return circuit
>>>
>>> # Under the unary encodings, occupation ``n`` of mode ``m`` sets qubit ``m * d + n``, with
>>> # ``d = max_occupation + 1`` levels per mode. That places |1, 1> on the basis state below.
>>> levels = max_occupation + 1
>>> index = sum(1 << (mode * levels + n) for mode, n in enumerate([1, 1]))
>>> initial = Statevector.from_int(index, 2**(num_sites * levels)).data
>>> exact = expm_multiply(-1j * total_time * matrix, initial)

>>> def error(circuit):
...     """Distance to the exact evolution, of the circuit's product of group exponentials."""
...     lowered = pm.run(circuit)
...     state = Statevector.from_int(0, 2**lowered.num_qubits).evolve(lowered).data
...     return float(np.max(np.abs(state - exact)))

Nota

error simula deliberadamente el circuito reducido sin descomponer sus instancias de PauliEvolutionGate. Statevector aplica cada una como la exponencial exacta de su grupo, de modo que lo único que queda es el producto de exponenciales de grupo que eligió la síntesis bosónica.

El valor predeterminado

Sin el argumento synthesis, una puerta Evolution se descompone con BosonicLieTrotter, la fórmula de primer orden: un factor por grupo, cada uno durante todo el tiempo.

>>> dict(build().decompose().count_ops())
{'Evolution': 2, 'InitializeModes': 1}

Sin descomponer, la puerta completa pasa a la etapa de qubits y no hay ningún error de Trotter bosónico:

>>> error(build()) < 1e-6
True

Al descomponer aparece un error real, porque los dos grupos no conmutan:

>>> error(build().decompose()) > 0.1
True

Reducir el error

Pasa un método de síntesis distinto a la puerta. Un orden mayor compone los factores de forma simétrica para que los términos principales del error se cancelen; más repeticiones dividen el tiempo en pasos más cortos.

>>> from qiskit_bosons.circuit.library.synthesis import BosonicLieTrotter, BosonicSuzukiTrotter
>>>
>>> for label, synthesis in [
...     ("Lie-Trotter", BosonicLieTrotter()),
...     ("Suzuki, order 2", BosonicSuzukiTrotter(order=2)),
...     ("Suzuki, order 4", BosonicSuzukiTrotter(order=4)),
...     ("Suzuki, order 4*", BosonicSuzukiTrotter(order=4, merge_adjacent=False)),
... ]:
...     circuit = build(synthesis=synthesis).decompose()
...     print(f"{label:16s} {circuit.count_ops()['Evolution']:3d} factors")
Lie-Trotter        2 factors
Suzuki, order 2    3 factors
Suzuki, order 4   11 factors
Suzuki, order 4*  15 factors

Cada aumento del orden en dos compone cinco barridos del orden inferior (el crecimiento $5^{(k-2)/2}$ de la recursión de Suzuki, aplicado a los $2n - 1$ factores del orden 2), que es el recuento de la última línea, marcado con una estrella. Los barridos coinciden en el mismo factor, de modo que las puertas vecinas evolucionan bajo un mismo operador y se fusionan en una sola puerta con el tiempo sumado. Eso ahorra una puerta por unión y es exacto. Pasa merge_adjacent=False para conservar la secuencia de los libros de texto. En cualquier caso, un orden mayor consigue un error menor:

>>> lie = error(build(synthesis=BosonicLieTrotter()).decompose())
>>> second = error(build(synthesis=BosonicSuzukiTrotter(order=2)).decompose())
>>> fourth = error(build(synthesis=BosonicSuzukiTrotter(order=4)).decompose())
>>> lie > second > fourth
True

Las repeticiones intercambian profundidad por precisión a un ritmo distinto. El error de primer orden decrece linealmente con el tamaño del paso, el de segundo orden cuadráticamente:

>>> lie_errors = [error(build(synthesis=BosonicLieTrotter(reps=r)).decompose()) for r in (1, 2, 4, 8)]
>>> bool(all(a > b for a, b in zip(lie_errors, lie_errors[1:])))
True
>>> bool(lie_errors[0] / lie_errors[-1] > 6)  # about 8, i.e. linear in 1 / reps
True
>>> suzuki_errors = [
...     error(build(synthesis=BosonicSuzukiTrotter(order=2, reps=r)).decompose()) for r in (1, 2, 4)
... ]
>>> bool(suzuki_errors[0] / suzuki_errors[-1] > 10)  # about 16, i.e. quadratic in 1 / reps
True

Cuál es el mejor compromiso depende del operador. Ten en cuenta que un operador con un solo grupo es exacto a cualquier orden, porque un factor aislado conmuta consigo mismo, y la fórmula se reduce entonces a una sola puerta en lugar de volver a dividirla.

Medir el coste

Un error menor cuesta más factores, y cada factor cuesta puertas de entrelazamiento una vez que la etapa de qubits lo ha sintetizado. Cuántas depende de la codificación, así que cuéntalas a través del pipeline preset, con la base fija y sin optimización a nivel de qubits (optimization_level=0), para ver el precio de cada método:

>>> basis = ["rz", "sx", "x", "cx"]
>>>
>>> def cx_count(circuit, encoding):
...     spec = LevelSpec(max_occupation, encoding)
...     lowered = generate_preset_occupation_pass_manager(
...         spec, optimization_level=0, basis_gates=basis
...     ).run(circuit)
...     return lowered.count_ops().get("cx", 0)
>>>
>>> methods = [
...     ("Lie-Trotter", BosonicLieTrotter()),
...     ("Suzuki, order 2", BosonicSuzukiTrotter(order=2)),
...     ("Suzuki, order 4", BosonicSuzukiTrotter(order=4)),
... ]
>>> encodings = ("gray", "unary_unprojected")
>>> counts = {
...     encoding: [
...         cx_count(build(synthesis=method).decompose(), encoding) for _, method in methods
...     ]
...     for encoding in encodings
... }
>>> # the order buys accuracy at a price that grows with it
>>> bool(counts["gray"][0] < counts["gray"][1] < counts["gray"][2])
True

Los recuentos exactos dependen de cómo sintetiza Qiskit las puertas de evolución, que cambia entre versiones, por lo que el gráfico los muestra en lugar de una tabla:

Cada grupo de barras es un método: Lie-Trotter, Suzuki de orden 2 y Suzuki de orden 4. El orden compra precisión a un precio que crece según el factor de la recursión, y la codificación unaria paga más que el código de Gray en todos los métodos porque cada uno de sus factores se asigna a más puertas de entrelazamiento. Elige el método según la precisión que necesites y la codificación según el hardware, y mantén separadas las dos decisiones: el circuito bosónico de arriba nunca cambia.

Aplicar un método a un circuito completo

Para establecer synthesis en cada puerta, tienes que transmitir la elección a través de todo lo que construye una. El pase B2BSynthesis toma la decisión una sola vez para un pipeline. Su config indica el método para cada tipo de puerta, y a la vez selecciona el método y expande las puertas:

>>> from qiskit_bosons.transpiler import BosonicPassManager
>>> from qiskit_bosons.transpiler.passes import B2BSynthesis
>>>
>>> trotterize = BosonicPassManager([B2BSynthesis({"Evolution": BosonicSuzukiTrotter(order=2)})])
>>> dict(trotterize.run(build()).count_ops())
{'Evolution': 3, 'InitializeModes': 1}

Un método también se puede indicar por nombre, con los argumentos de su __init__, en lugar de construirlo: B2BSynthesis({"Evolution": ("SuzukiTrotter", (), {"order": 2})}). Los nombres disponibles son los de B2BSynthesisPluginManager, y un método que no está registrado en ningún sitio se puede entregar como instancia, como arriba.

Con apply=False el pase solo selecciona el método y deja la expansión para el momento en que se quiera, por ejemplo para inspeccionar antes el circuito. Una puerta atómica nunca se toca, ni siquiera se etiqueta:

>>> selected = BosonicPassManager(
...     [B2BSynthesis({"Evolution": ("SuzukiTrotter", (), {"order": 2})}, apply=False)]
... ).run(build())
>>> selected.count_ops()["Evolution"]  # selected, not yet expanded
1
>>> selected.decompose().count_ops()["Evolution"]  # expanded on demand
3

Una puerta sin expandir llega entera a la etapa de bosón a qubit, donde se mapea sin consultar su synthesis, de modo que todos los métodos producen entonces una salida idéntica.

Un filter opcional lo restringe a las puertas que se benefician. Las evoluciones cuyos grupos conmutan todos, como un hamiltoniano formado solo por operadores de número, se sintetizan exactamente a cualquier orden, así que un orden mayor solo añadiría profundidad ahí.

Restringir el pase a las puertas que se benefician

Un operador puramente diagonal hace concreto el argumento, porque los operadores de número conmutan entre sí. Sus grupos son exactos en primer orden, y el segundo orden gasta puertas de entrelazamiento para nada:

>>> # one group per term, each Hermitian; the order is fixed because order two evolves the first
>>> # group twice, so which group comes first changes the cost
>>> diagonal = BosonOperator.from_terms_with_groups(
...     [([cre(0), ann(0), cre(1), ann(1)], 1.0, 0), ([cre(0), ann(0)], 0.5, 1)]
... )
>>>
>>> def diagonal_cost(synthesis):
...     circuit = BosonicCircuit(num_sites)
...     circuit.append(Evolution(num_sites, diagonal, time=total_time, synthesis=synthesis), circuit.modes)
...     return cx_count(circuit.decompose(), "gray")
>>>
>>> bool(diagonal_cost(BosonicSuzukiTrotter(order=2)) > diagonal_cost(BosonicLieTrotter()))
True

Ambos circuitos son exactos, así que el segundo gasta más puertas de entrelazamiento sin ganar nada. (La fórmula simétrica de orden dos aplica el primer grupo dos veces con la mitad del tiempo cada una, y por eso el coste depende de qué grupo va primero.) El filter toma un predicado sobre DAGOpNode y deja intactas las puertas rechazadas, de modo que un solo pase puede subir el orden en las puertas que se benefician y mantener el resto en primer orden. El predicado siguiente acepta una puerta solo cuando su operador mueve bosones entre modos, que es cuando los grupos no conmutan:

>>> def is_diagonal(operator):
...     """Whether every term only counts bosons, 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 = BosonicCircuit(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
>>>
>>> order_2 = BosonicSuzukiTrotter(order=2)
>>> everywhere = cx_count(
...     BosonicPassManager([B2BSynthesis({"Evolution": order_2})]).run(mixed_circuit()), "gray"
... )
>>> restricted = cx_count(
...     BosonicPassManager(
...         [
...             B2BSynthesis(
...                 {"Evolution": order_2},
...                 filter=lambda node: not is_diagonal(node.op.operator),
...             )
...         ]
...     ).run(mixed_circuit()),
...     "gray",
... )
>>> bool(everywhere > restricted)  # the filter saves entangling gates
True

Si vuelves a ejecutar el pase, no cambia nada: los factores que emite son atómicos.

Factores atómicos

Una Evolution marcada como atomic no tiene definición, así que decompose() la deja donde está por muchas veces que se aplique. Una síntesis marca así los factores que emite, y también lo hace BosonicQDrift con los términos que muestrea.

Esto es más que orden. Dividir un grupo hermítico más allá, término a término, produce factores que individualmente no son hermíticos, y la exponencial de un operador no hermítico no es unitaria:

>>> circuit = build()
>>> once = circuit.decompose()
>>> dict(circuit.decompose(reps=5).count_ops()) == dict(once.count_ops())
True

Precaución

Cada factor debe ser hermítico por sí mismo. Un operador hermítico dividido en grupos puede producir grupos no hermíticos, como un término de salto separado de su adjunto. Compruébalo con groups_are_hermitian(). Un operador sin grupos se divide término a término, y la etapa de qubits rechaza entonces los factores con Operator contains complex coefficients. La simetrización de un orden par cancela parcialmente el error de un factor no hermítico, de modo que un operador mal agrupado puede parecer mejor en el orden 2 que en el orden 1 y seguir siendo incorrecto.

Repetir un paso corto

Como alternativa a nivel de circuito, construye un paso de Trotter corto y repítelo. repeat() aplana las repeticiones en un solo circuito y deja los límites adyacentes, de modo que los pases pueden actuar a través de ellos:

>>> step = BosonicCircuit(num_sites)
>>> step.append(Evolution(num_sites, hamiltonian, time=total_time / 4), step.modes)
>>> dict(step.repeat(4).count_ops())
{'Evolution': 4}

Pasa insert_barriers=True para mantener separadas las repeticiones. Una barrera se mapea a todos los qubits de los registros que toca, porque una codificación necesita varios qubits por modo y un modo no tiene una imagen propia por qubit:

>>> [inst.operation.name for inst in step.repeat(3, insert_barriers=True)._inner.data]
['Evolution', 'barrier', 'Evolution', 'barrier', 'Evolution']

Nota

repeat no guarda relación con reps de un método de síntesis. Una síntesis subdivide una única evolución en más factores, más cortos. Repetir un circuito que contiene una evolución simplemente multiplica su tiempo total.

Combinar con una fórmula a nivel de qubits

Los dos niveles se componen, y donde ambos aproximan, la precisión del resultado la rige el más débil de los dos. Subir el orden bosónico mientras la fórmula del lado de los qubits sigue en primer orden aporta poco. Consulta MapperFnEvolutionSynthesis para la fórmula del lado de los qubits. Conserva el orden de los términos de Pauli que eligió el mapeador, y así group_wise() puede emitir los términos de cada grupo adyacentes entre sí, y simplify() puede fijar en su lugar un orden canónico.

Siguientes pasos

• Transpilar circuitos bosónicos explica las etapas que siguen a la síntesis bosónica, y dónde se sitúa B2BSynthesis entre ellas.

• Agrupar términos del operador: usar la estructura del operador explica cómo asignar y comprobar los grupos que entrelaza una fórmula de producto.

• Simular la dinámica de Bose-Hubbard con una codificación a elegir aplica el esquema de repetición a un ejemplo completo.