Synthétiser une porte d'évolution dans l'espace bosonique
L'évolution temporelle d'un hamiltonien qui est une somme de morceaux non commutants ne peut pas être implémentée exactement porte par porte. Une formule de produit l'approxime par une suite de courtes évolutions sous chacun des morceaux, et l'erreur d'approximation s'appelle l'erreur de Trotter.
Tu peux appliquer une formule de produit à deux endroits, et ils ne sont pas interchangeables :
-
Dans l'espace bosonique, avant qu'une troncature ou un encodage n'ait été choisi. La formule voit les groupes de l'opérateur (voir Grouper les termes d'un opérateur : exploiter la structure de l'opérateur) et peut les entrelacer.
-
Dans l'espace des qubits, à l'intérieur de la
PauliEvolutionGatede Qiskit, une fois l'opérateur mappé. La formule voit un opérateur de Pauli par groupe et ne peut subdiviser que les termes à l'intérieur de celui-ci.
Seule la première peut réduire l'erreur entre des groupes qui ne commutent pas, qui est généralement l'erreur dominante. Ce guide montre comment l'utiliser.
Les deux étapes sont facultatives. Une porte Evolution qui n'est jamais décomposée dans l'espace bosonique atteint entière l'étape boson-vers-qubit, où elle est mappée et synthétisée. Cela n'entraîne aucune erreur au niveau bosonique, seulement celle qu'introduit la formule au niveau des qubits.
Configuration
Une chaîne de Bose-Hubbard à deux sites, le plus petit exemple où les groupes ne commutent pas (voir Simuler la dynamique de Bose-Hubbard avec un choix d'encodage pour le modèle complet). Elle est groupée en un groupe hermitien pour la liaison et un groupe diagonal pour l'interaction sur 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
Pour mesurer une erreur, il faut une référence. L'évolution exacte du hamiltonien tel que mappé à la même troncature est la bonne, car la troncature est une approximation distincte :
>>> 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)))
error simule délibérément le circuit abaissé sans décomposer ses instances de PauliEvolutionGate. Statevector applique chacune comme l'exponentielle exacte de son groupe, de sorte que la seule chose qui reste est le produit d'exponentielles de groupes choisi par la synthèse bosonique.
Le défaut
Sans argument synthesis, une porte Evolution se décompose avec BosonicLieTrotter, la formule du premier ordre : un facteur par groupe, chacun pour la durée totale.
>>> dict(build().decompose().count_ops())
{'Evolution': 2, 'InitializeModes': 1}
Sans décomposition, la porte entière est transmise à l'étape qubit et il n'y a aucune erreur de Trotter bosonique :
>>> error(build()) < 1e-6
True
Une fois décomposée, une véritable erreur apparaît, car les deux groupes ne commutent pas :
>>> error(build().decompose()) > 0.1
True
Réduire l'erreur
Passe une autre méthode de synthèse à la porte. Un ordre plus élevé compose les facteurs de façon symétrique pour que les termes d'erreur dominants s'annulent ; davantage de répétitions divisent le temps en pas plus courts.
>>> 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
Chaque augmentation de l'ordre de deux compose cinq balayages de l'ordre inférieur (la croissance en de la récursion de Suzuki, appliquée aux facteurs de l'ordre 2), ce qui correspond au décompte de la dernière ligne, marqué d'une étoile. Les balayages se rejoignent sur le même facteur, de sorte que les portes voisines évoluent sous un même opérateur et sont fusionnées en une seule porte dont le temps est la somme. Cela économise une porte par jonction et reste exact. Passe merge_adjacent=False pour conserver la séquence théorique. Dans les deux cas, un ordre plus élevé permet d'obtenir une erreur plus petite :
>>> 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
Les répétitions échangent la profondeur contre la précision à un rythme différent. L'erreur du premier ordre décroît linéairement avec la taille du pas, celle du second ordre quadratiquement :
>>> 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
Le meilleur compromis dépend de l'opérateur. Note qu'un opérateur à un seul groupe est exact à n'importe quel ordre, car un facteur isolé commute avec lui-même, et la formule se réduit alors à une seule porte au lieu de la subdiviser de nouveau.
Mesurer le coût
Une erreur plus faible coûte plus de facteurs, et chaque facteur coûte des portes intriquantes une fois synthétisé par l'étape qubit. Leur nombre dépend de l'encodage ; compte-les donc via le pipeline du preset, avec la base fixée et sans optimisation au niveau des qubits (optimization_level=0), pour voir le prix de chaque méthode :
>>> 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
Les décomptes exacts dépendent de la façon dont Qiskit synthétise les portes d'évolution, ce qui change d'une version à l'autre ; le graphique les montre donc plutôt qu'un tableau :
Chaque groupe de barres correspond à une méthode : Lie-Trotter, Suzuki d'ordre 2 et Suzuki d'ordre 4. L'ordre achète de la précision à un prix qui croît du facteur de la récursion, et l'encodage unaire paie plus que le code de Gray pour chaque méthode parce que chacun de ses facteurs se traduit par plus de portes intriquantes. Choisis la méthode selon la précision dont tu as besoin et l'encodage selon le matériel, et garde les deux choix séparés : le circuit bosonique ci-dessus ne change jamais.
Appliquer une méthode à un circuit entier
Pour définir synthesis sur chaque porte, tu dois faire passer ce choix à travers tout ce qui en construit une. La passe B2BSynthesis fait ce choix une seule fois pour un pipeline. Son config nomme la méthode pour chaque type de porte, et elle à la fois sélectionne la méthode et développe les portes :
>>> 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}
Une méthode peut aussi être désignée par son nom, avec les arguments de son __init__, au lieu d'être construite : B2BSynthesis({"Evolution": ("SuzukiTrotter", (), {"order": 2})}). Les noms disponibles sont ceux du B2BSynthesisPluginManager, et une méthode enregistrée nulle part peut être transmise comme instance, comme ci-dessus.
Avec apply=False, la passe se contente de sélectionner la méthode et laisse le développement au moment où il est voulu, par exemple pour inspecter d'abord le circuit. Une porte atomique n'est jamais touchée, pas même étiquetée :
>>> 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
Une porte non développée atteint entière l'étape boson-vers-qubit, où elle est mappée sans que son synthesis soit consulté ; chaque méthode produit alors un résultat identique.
Un filter facultatif la restreint aux portes qui en profitent. Les évolutions dont tous les groupes commutent, comme un hamiltonien composé uniquement d'opérateurs nombre, sont synthétisées exactement à n'importe quel ordre ; un ordre plus élevé n'y ajouterait donc que de la profondeur.
Restreindre la passe aux portes qui en profitent
Un opérateur purement diagonal rend le propos concret, car les opérateurs nombre commutent entre eux. Ses groupes sont exacts au premier ordre, et le second ordre dépense des portes intriquantes pour rien :
>>> # 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
Les deux circuits sont exacts ; le second dépense donc plus de portes intriquantes sans rien gagner. (La formule symétrique d'ordre deux applique le premier groupe deux fois à la moitié du temps chacune, et le coût dépend donc du groupe qui vient en premier.) Le filter prend un prédicat sur le DAGOpNode et laisse les portes rejetées intactes, de sorte qu'une même passe peut élever l'ordre sur les portes qui en profitent et garder le reste au premier ordre. Le prédicat ci-dessous n'accepte une porte que si son opérateur déplace des bosons entre modes, ce qui est le cas lorsque les groupes ne commutent pas :
>>> 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 tu exécutes à nouveau la passe, rien ne change : les facteurs qu'elle émet sont atomiques.
Facteurs atomiques
Une Evolution marquée atomic n'a pas de définition, donc decompose() la laisse en place, quel que soit le nombre d'applications. Une synthèse marque ainsi les facteurs qu'elle émet, et BosonicQDrift fait de même pour les termes qu'il échantillonne.
C'est plus qu'une question d'ordre. Subdiviser davantage un groupe hermitien, terme par terme, donne des facteurs qui ne sont individuellement pas hermitiens, et l'exponentielle d'un opérateur non hermitien n'est pas unitaire :
>>> circuit = build()
>>> once = circuit.decompose()
>>> dict(circuit.decompose(reps=5).count_ops()) == dict(once.count_ops())
True
Chaque facteur doit être hermitien par lui-même. Un opérateur hermitien découpé en groupes peut produire des groupes non hermitiens, comme un terme de saut séparé de son adjoint. Vérifie avec groups_are_hermitian(). Un opérateur sans groupes est découpé terme par terme, et l'étape qubit rejette alors les facteurs avec Operator contains complex coefficients. La symétrisation d'un ordre pair annule partiellement l'erreur d'un facteur non hermitien ; un opérateur mal groupé peut donc sembler meilleur à l'ordre 2 qu'à l'ordre 1 tout en restant faux.
Répéter un pas court
Pour une alternative au niveau du circuit, construis un court pas de Trotter et répète-le. repeat() aplatit les répétitions en un seul circuit et laisse les frontières adjacentes, de sorte que les passes peuvent agir à travers elles :
>>> step = BosonicCircuit(num_sites)
>>> step.append(Evolution(num_sites, hamiltonian, time=total_time / 4), step.modes)
>>> dict(step.repeat(4).count_ops())
{'Evolution': 4}
Passe insert_barriers=True pour garder les répétitions séparées. Une barrière est mappée sur chaque qubit des registres qu'elle touche, parce qu'un encodage a besoin de plusieurs qubits par mode et qu'un mode n'a pas d'image par qubit qui lui soit propre :
>>> [inst.operation.name for inst in step.repeat(3, insert_barriers=True)._inner.data]
['Evolution', 'barrier', 'Evolution', 'barrier', 'Evolution']
repeat n'a aucun rapport avec les reps d'une méthode de synthèse. Une synthèse subdivise une seule évolution en facteurs plus nombreux et plus courts. Répéter un circuit qui contient une évolution ne fait que multiplier son temps total.
Combiner avec une formule au niveau des qubits
Les deux niveaux se composent, et là où les deux approximent, la précision du résultat est dictée par le plus faible des deux. Élever l'ordre bosonique alors que la formule côté qubit reste au premier ordre apporte peu. Voir MapperFnEvolutionSynthesis pour la formule côté qubit. Elle conserve l'ordre des termes de Pauli choisi par le mapper, et ainsi group_wise() peut émettre les termes de chaque groupe adjacents les uns aux autres, tandis que simplify() peut imposer à la place un ordre canonique.
Étapes suivantes
-
Transpiler des circuits bosoniques explique les étapes qui suivent la synthèse bosonique, et où
B2BSynthesisse situe parmi elles. -
Grouper les termes d'un opérateur : exploiter la structure de l'opérateur explique comment attribuer et vérifier les groupes qu'une formule de produit entrelace.
-
Simuler la dynamique de Bose-Hubbard avec un choix d'encodage applique le schéma de répétition à un exemple complet.