Sintetizzare un gate Evolution nello spazio bosonico

L'evoluzione temporale $e^{-itH}$ di una Hamiltoniana che è una somma di pezzi non commutanti non può essere implementata esattamente gate per gate. Una formula di prodotto la approssima con una sequenza di brevi evoluzioni sotto i singoli pezzi, e l'errore di approssimazione si chiama errore di Trotter.

Puoi applicare una formula di prodotto in due punti, e non sono intercambiabili:

1. Nello spazio bosonico, prima che sia stato scelto qualsiasi troncamento o codifica. La formula vede i gruppi dell'operatore (vedi Raggruppare i termini dell'operatore: usare la struttura dell'operatore) e può intercalarli.

2. Nello spazio dei qubit, all'interno del PauliEvolutionGate di Qiskit, dopo che l'operatore è stato mappato. La formula vede un operatore di Pauli per gruppo e può solo suddividere i termini al suo interno.

Solo la prima può ridurre l'errore tra gruppi che non commutano, che di solito è quello dominante. Questa guida mostra come usarla.

Importante

Entrambi i passaggi sono facoltativi. Un gate Evolution che non viene mai decomposto nello spazio bosonico arriva intero alla fase da bosone a qubit e viene mappato e sintetizzato lì. Ciò non comporta alcun errore a livello bosonico, solo quello che introduce la formula a livello di qubit.

Configurazione

Una catena di Bose-Hubbard a due siti, l'esempio più piccolo in cui i gruppi non commutano (vedi Simulare la dinamica di Bose-Hubbard con una scelta di codifica per il modello completo). È raggruppata in un gruppo hermitiano per il legame e un gruppo diagonale per l'interazione on-site:

>>> 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

Per misurare un errore serve un riferimento. Quello giusto è l'evoluzione esatta dell'Hamiltoniana così come mappata con lo stesso troncamento, perché il troncamento è un'approssimazione separata:

>>> 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 deliberatamente il circuito tradotto senza decomporre le sue istanze di PauliEvolutionGate. Statevector applica ciascuna come l'esponenziale esatto del suo gruppo, quindi l'unica cosa che resta è il prodotto di esponenziali di gruppo scelto dalla sintesi bosonica.

Il comportamento predefinito

Senza l'argomento synthesis, un gate Evolution si decompone con BosonicLieTrotter, la formula del primo ordine: un fattore per gruppo, ciascuno per l'intero tempo.

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

Senza decomposizione, l'intero gate viene passato alla fase dei qubit e non c'è alcun errore di Trotter bosonico.

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

Dopo la decomposizione compare un errore reale, perché i due gruppi non commutano:

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

Ridurre l'errore

Passa al gate un metodo di sintesi diverso. Un ordine più alto compone i fattori in modo simmetrico così che i termini d'errore dominanti si cancellino; più ripetizioni dividono il tempo in passi più brevi.

>>> 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

Ogni aumento dell'ordine di due compone cinque passate di quello sottostante (la crescita $5^{(k-2)/2}$ della ricorsione di Suzuki, applicata ai $2n - 1$ fattori dell'ordine 2), che è il conteggio nell'ultima riga, contrassegnato da una stella. Le passate si incontrano sullo stesso fattore, quindi i gate adiacenti evolvono sotto un solo operatore e vengono fusi in un unico gate con il tempo sommato. Ciò fa risparmiare un gate per giuntura ed è esatto. Passa merge_adjacent=False per mantenere la sequenza da manuale. In entrambi i casi, un ordine più alto fa ottenere un errore minore:

>>> 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

Le ripetizioni scambiano profondità con accuratezza a un tasso diverso. L'errore del primo ordine diminuisce linearmente con la dimensione del passo, quello del secondo ordine quadraticamente:

>>> 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

Quale sia lo scambio migliore dipende dall'operatore. Nota che un operatore con un solo gruppo è esatto a qualsiasi ordine, perché un fattore isolato commuta con se stesso, e la formula si riduce quindi a un solo gate invece di suddividerlo di nuovo.

Misurare il costo

Un errore minore costa più fattori, e ogni fattore costa gate di entanglement una volta che la fase dei qubit lo ha sintetizzato. Quanti dipende dalla codifica, quindi contali attraverso la pipeline del preset, con la base fissata e senza ottimizzazione a livello di qubit (optimization_level=0), per vedere il prezzo di ciascun metodo:

>>> 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

I conteggi esatti dipendono da come Qiskit sintetizza i gate di evoluzione, che cambia tra una release e l'altra, quindi il grafico li mostra al posto di una tabella:

Ogni gruppo di barre è un metodo: Lie-Trotter, Suzuki di ordine 2 e Suzuki di ordine 4. L'ordine compra accuratezza a un prezzo che cresce con il fattore della ricorsione, e la codifica unaria paga più del codice di Gray per ogni metodo perché ciascuno dei suoi fattori si mappa su più gate di entanglement. Scegli il metodo per l'accuratezza che ti serve e la codifica per l'hardware, e tieni separate le due scelte: il circuito bosonico qui sopra non cambia mai.

Applicare un metodo a un intero circuito

Per impostare synthesis su ogni gate, devi propagare la scelta attraverso tutto ciò che ne costruisce uno. Il pass B2BSynthesis compie invece la scelta una sola volta per una pipeline. Il suo config indica il metodo per ogni tipo di gate, e sia seleziona il metodo sia espande i gate:

>>> 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 metodo può anche essere indicato per nome, con gli argomenti per il suo __init__, invece di essere costruito: B2BSynthesis({"Evolution": ("SuzukiTrotter", (), {"order": 2})}). I nomi disponibili sono quelli del B2BSynthesisPluginManager, e un metodo non registrato da nessuna parte può essere passato come istanza, come sopra.

Con apply=False il pass seleziona soltanto il metodo e lascia l'espansione al momento in cui serve, per esempio per ispezionare prima il circuito. Un gate atomico non viene mai toccato, nemmeno etichettato:

>>> 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

Un gate non espanso arriva intero alla fase da bosone a qubit, dove viene mappato senza che il suo synthesis venga consultato, quindi ogni metodo produce allora un output identico.

Un filter facoltativo lo limita ai gate che ne traggono vantaggio. Le evoluzioni i cui gruppi commutano tutti, come una Hamiltoniana fatta solo di operatori di numero, vengono sintetizzate esattamente a qualsiasi ordine, quindi un ordine più alto vi aggiungerebbe soltanto profondità.

Limitare il pass ai gate che ne traggono vantaggio

Un operatore puramente diagonale rende il punto concreto, perché gli operatori di numero commutano tra loro. I suoi gruppi sono esatti al primo ordine, e il secondo ordine spende gate di entanglement per nulla:

>>> # 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

Entrambi i circuiti sono esatti, quindi il secondo spende più gate di entanglement senza alcun guadagno. (La formula simmetrica di ordine due applica il primo gruppo due volte a metà tempo ciascuna, e quindi il costo dipende da quale gruppo viene per primo.) Il filter prende un predicato sul DAGOpNode e lascia intatti i gate rifiutati, così un solo pass può alzare l'ordine sui gate che ne traggono vantaggio e mantenere gli altri al primo ordine. Il predicato qui sotto accetta un gate solo quando il suo operatore sposta bosoni tra modi, che è il caso in cui i gruppi non commutano:

>>> 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

Se esegui di nuovo il pass, non cambia nulla: i fattori che emette sono atomici.

Fattori atomici

Un Evolution contrassegnato come atomic non ha definizione, quindi decompose() lo lascia al suo posto per quante volte venga applicato. Una sintesi contrassegna in questo modo i fattori che emette, e lo fa anche BosonicQDrift per i termini che campiona.

Questo è più che una questione di ordine. Suddividere ulteriormente un gruppo hermitiano, termine per termine, produce fattori che individualmente non sono hermitiani, e l'esponenziale di un operatore non hermitiano non è unitario:

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

Attenzione

Ogni fattore deve essere hermitiano da solo. Un operatore hermitiano suddiviso in gruppi può produrne di non hermitiani, come un termine di hopping separato dal suo aggiunto. Verifica con groups_are_hermitian(). Un operatore senza gruppi viene suddiviso termine per termine, e la fase dei qubit rifiuta poi i fattori con Operator contains complex coefficients. La simmetrizzazione di un ordine pari cancella in parte l'errore di un fattore non hermitiano, quindi un operatore raggruppato in modo errato può sembrare migliore all'ordine 2 che all'ordine 1 pur essendo sbagliato.

Ripetere un passo breve

Per un'alternativa a livello di circuito, costruisci un breve passo di Trotter e ripetilo. repeat() appiattisce le ripetizioni in un unico circuito e lascia adiacenti i confini, così i pass possono agire attraverso di essi:

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

Passa insert_barriers=True per mantenere separate le ripetizioni. Una barriera viene mappata su ogni qubit dei registri che tocca, perché una codifica richiede più qubit per modo e un modo non ha un'immagine per singolo qubit propria:

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

Nota

repeat non ha relazione con i reps di un metodo di sintesi. Una sintesi suddivide una singola evoluzione in più fattori, più brevi. Ripetere un circuito che contiene una sola evoluzione si limita a moltiplicarne il tempo totale.

Combinare con una formula a livello di qubit

I due livelli si compongono, e dove entrambi approssimano, l'accuratezza del risultato è determinata dal più debole dei due. Alzare l'ordine bosonico mentre la formula lato qubit resta al primo ordine serve a poco. Vedi MapperFnEvolutionSynthesis per la formula lato qubit. Mantiene l'ordine dei termini di Pauli scelto dal mapper, e quindi group_wise() può emettere i termini di ciascun gruppo adiacenti tra loro, e simplify() può invece fissare un ordine canonico.

Prossimi passi

• Eseguire il transpile di circuiti bosonici spiega le fasi che seguono la sintesi bosonica, e dove si colloca tra esse B2BSynthesis.

• Raggruppare i termini dell'operatore: usare la struttura dell'operatore spiega come assegnare e verificare i gruppi che una formula di prodotto intercala.

• Simulare la dinamica di Bose-Hubbard con una scelta di codifica applica lo schema di ripetizione a un esempio completo.