Générer des circuits d'évolution temporelle de Krylov (SKQD)

Important

Les concepts de ce guide ne sont actuellement disponibles que dans l'API Python. Une fonctionnalité équivalente sera disponible dans l'API C dans une prochaine version.

La diagonalisation quantique de Krylov basée sur l'échantillonnage (SKQD) est une variante quantique-centrée de la diagonalisation en sous-espace. Plutôt que d'échantillonner un unique ansatz variationnel (comme dans la diagonalisation quantique basée sur l'échantillonnage (SQD) simple) ou un ensemble de circuits d'évolution temporelle aléatoires (comme dans SqDRIFT), elle échantillonne une base de Krylov : la famille d'états

$$\lvert \psi_k \rangle = e^{-i k \, \Delta t \, H} \lvert \psi_0 \rangle, \qquad k = 0, 1, \dots, D - 1,$$

obtenue en faisant évoluer un état de référence $\lvert \psi_0 \rangle$ pendant des multiples croissants d'un pas de temps fixe $\Delta t$. Les chaînes de bits échantillonnées à partir de ces $D$ circuits peuplent le sous-espace de Krylov dans lequel l'hamiltonien est ensuite diagonalisé classiquement.

Ce guide reproduit l'étape de construction des circuits de l'algorithme SKQD pour le modèle d'Anderson à impureté unique (SIAM), une impureté magnétique couplée à un bain sans interaction, en suivant la publication SKQD.

Remarque

L'étape d'échantillonnage de SKQD ne fonctionne que si l'état fondamental est creux dans la base computationnelle : la diagonalisation classique doit pouvoir le reconstruire à partir d'un nombre raisonnable de chaînes de bits échantillonnées. Comme on le voit plus bas, cette parcimonie n'est pas automatique ; c'est une propriété de la base à particule unique dans laquelle le problème est exprimé. Bien choisir la base est le point crucial de cette construction.

1. Construire l'hamiltonien SIAM comme opérateur fermionique

Le modèle d'Anderson à impureté unique place une orbitale d'impureté en interaction (indice $0$) en contact avec un bain sans interaction de $L$ sites disposés en chaîne. Son hamiltonien se décompose en une partie à un corps $H_1$ et une partie à deux corps $H_2$,

$$H = H_1 + H_2,$$

la partie à un corps rassemblant l'énergie locale de l'impureté (le potentiel chimique $\mu$), l'hybridation impureté-bain $V$ et le saut du bain entre plus proches voisins $t$,

$$H_1 = \mu \sum_\sigma a^\dagger_{0\sigma} a_{0\sigma} - V \sum_\sigma \left( a^\dagger_{0\sigma} a_{1\sigma} + \text{h.c.} \right) - t \sum_{\sigma} \sum_{i=1}^{L-1} \left( a^\dagger_{i\sigma} a_{i+1,\sigma} + \text{h.c.} \right),$$

et la partie à deux corps portant la répulsion de Coulomb locale $U$ sur la seule impureté,

$$H_2 = U \, a^\dagger_{0\uparrow} a_{0\uparrow} a^\dagger_{0\downarrow} a_{0\downarrow}.$$

Ce guide se place au point de symétrie particule-trou $\mu = -U/2$. Ces deux parties correspondent directement aux intégrales à un corps h1e (saut du bain, hybridation de l'impureté, énergie locale de l'impureté) et à l'unique intégrale à deux corps h2e (le terme de Coulomb de l'impureté). Assemble-les d'abord dans la base de position, où le modèle est naturellement défini.

>>> import numpy as np
>>>
>>> norb = 8
>>> nelec = (norb // 2, norb // 2)  # half filling
>>> hopping, onsite, hybridization = 1.0, 10.0, 1.0
>>> chemical_potential = -0.5 * onsite
>>>
>>> # one-body integrals: bath hopping, impurity hybridization, and impurity on-site energy
>>> h1e = np.zeros((norb, norb))
>>> np.fill_diagonal(h1e[:, 1:], -hopping)
>>> np.fill_diagonal(h1e[1:, :], -hopping)
>>> h1e[0, 1] = h1e[1, 0] = -hybridization
>>> h1e[0, 0] = chemical_potential
>>>
>>> # two-body integrals: on-site Coulomb repulsion on the impurity orbital
>>> h2e = np.zeros((norb,) * 4)
>>> h2e[0, 0, 0, 0] = onsite

2. Passer à la base des impulsions

La partie bain du SIAM est une chaîne de fermions libres, elle est donc diagonalisée par un changement vers la base des impulsions. C'est l'étape décisive pour SKQD. Comme le bain est invariant par translation, l'état fondamental n'est pas creux dans la base de position, mais il l'est dans la base des impulsions, où l'hamiltonien du bain est diagonal. La rotation d'orbitales à particule unique qui effectue ce changement relève de l'algèbre linéaire propre au problème. Elle diagonalise la matrice de saut du bain en laissant l'orbitale d'impureté intacte, puis déplace l'impureté vers un site central :

>>> def momentum_basis(norb):
...     """Orbital rotation from the position to the momentum basis."""
...     n_bath = norb - 1
...     hopping_matrix = np.zeros((n_bath, n_bath))
...     np.fill_diagonal(hopping_matrix[:, 1:], -1)
...     np.fill_diagonal(hopping_matrix[1:, :], -1)
...     _, vecs = np.linalg.eigh(hopping_matrix)
...     orbital_rotation = np.zeros((norb, norb))
...     orbital_rotation[0, 0] = 1.0
...     orbital_rotation[1:, 1:] = vecs
...     new_index = n_bath // 2
...     perm = np.r_[1:(new_index + 1), 0, (new_index + 1):norb]
...     return orbital_rotation[:, perm]
...
>>> orbital_rotation = momentum_basis(norb)

Contrairement à un workflow variationnel, on fait ici tourner les intégrales de l'hamiltonien vers la base des impulsions, de sorte que tout le circuit (hamiltonien comme état) vit dans la base où l'état fondamental est creux et où les chaînes de bits échantillonnées sont donc informatives. La rotation des intégrales à un et deux corps est un changement de base standard des tenseurs de structure électronique :

>>> rotation = orbital_rotation.T.conj()
>>> h1e_momentum = np.einsum("ab,Aa,Bb->AB", h1e, rotation, rotation.conj(), optimize="greedy")
>>> h2e_momentum = np.einsum(
...     "abcd,Aa,Bb,Cc,Dd->ABCD",
...     h2e,
...     rotation,
...     rotation.conj(),
...     rotation,
...     rotation.conj(),
...     optimize="greedy",
... )

Comme la rotation laisse l'orbitale d'impureté fixe (elle n'agit que sur le bain), l'interaction locale est inchangée. h2e_momentum reste un unique terme nombre-nombre, tout comme h2e. Le seul changement est que l'impureté a été déplacée vers un indice central. C'est visible directement dans les tenseurs d'intégrales, qui portent chacun une seule entrée non nulle :

>>> print(f"h2e nonzero at:          {np.argwhere(np.abs(h2e) > 1e-12).tolist()}")
h2e nonzero at:          [[0, 0, 0, 0]]
>>> print(f"h2e_momentum nonzero at: {np.argwhere(np.abs(h2e_momentum) > 1e-12).tolist()}")
h2e_momentum nonzero at: [[3, 3, 3, 3]]

Cette localité est ce qui garde l'interaction peu coûteuse à synthétiser par la suite.

Assemble les intégrales à un corps en un FermionOperator avec le constructeur d'intégrales électroniques from_1body_tril_spin_sym(), qui les attend empaquetées selon l'ordre triangulaire inférieur. La partie à deux corps est un unique terme nombre-nombre ; plutôt que de la faire passer par la mécanique des intégrales à deux corps, construis-la directement avec cre()/ann() : l'opérateur de Coulomb local $U \, n_{p\uparrow} n_{p\downarrow}$ pour le mode d'impureté $p$ (son mode $\uparrow$ à l'indice p, son mode $\downarrow$ à l'indice p + norb dans la disposition par blocs de spin). Les deux constructeurs doublent le spin des norb orbitales spatiales en 2 * norb modes fermioniques.

Garde la partie à un corps $H_1$ et la partie à deux corps $H_2$ comme des opérateurs séparés plutôt que de les additionner immédiatement, conformément aux termes $H_1$ et $H_2$ des équations de la section 1. Construire chaque partie une seule fois (dans les deux bases) permet de les réutiliser plus bas. Additionnées, elles donnent l'hamiltonien complet pour la diagonalisation exacte dans l'une ou l'autre base. Prise séparément, $H_2$ (base des impulsions) est l'opérateur que le pas de Trotter fait évoluer à l'étape 4.

>>> from qiskit_fermions.operators import FermionOperator, ann, cre
>>>
>>> # U n_imp_up n_imp_down: the impurity's up/down modes sit at `imp` and `imp + norb` in the
>>> # block-spin layout (the impurity is orbital 0 in the position basis, orbital 3 in the
>>> # momentum basis -- the single relocated index found above)
>>> imp_pos, imp_mom = 0, 3
>>>
>>> one_body_position = FermionOperator.from_1body_tril_spin_sym(
...     h1e[np.tril_indices(norb)], norb
... )
>>> two_body_position = FermionOperator.from_dict(
...     {(cre(imp_pos), ann(imp_pos), cre(imp_pos + norb), ann(imp_pos + norb)): onsite}
... )
>>> one_body_momentum = FermionOperator.from_1body_tril_spin_sym(
...     h1e_momentum[np.tril_indices(norb)], norb
... )
>>> two_body_momentum = FermionOperator.from_dict(
...     {(cre(imp_mom), ann(imp_mom), cre(imp_mom + norb), ann(imp_mom + norb)): onsite}
... )
>>>
>>> hamiltonian_position = (
...     (one_body_position + two_body_position).normal_ordered().simplify(atol=1e-14)
... )
>>> hamiltonian_momentum = (
...     (one_body_momentum + two_body_momentum).normal_ordered().simplify(atol=1e-14)
... )

Le changement de base est une transformation de similitude unitaire, il laisse donc le spectre inchangé. Confirme que l'opérateur est correct en comparant sa plus basse valeur propre dans le secteur $(\text{norb}, \text{nelec})$ à une diagonalisation exacte. Le FermionOperator expose un LinearOperator de SciPy (reposant sur un noyau natif de produit matrice-vecteur d'interaction de configurations complète (FCI)), il peut donc être passé directement à scipy.sparse.linalg.eigsh().

>>> import ffsim
>>> import scipy.sparse.linalg
>>>
>>> linop = ffsim.linear_operator(hamiltonian_momentum, norb, nelec)
>>> reference_energy, ground_state = scipy.sparse.linalg.eigsh(linop, k=1, which="SA")
>>> reference_energy = reference_energy[0]
>>> ground_state = ground_state[:, 0]
>>> print(f"exact ground-state energy: {reference_energy:.6f}")
exact ground-state energy: -13.422492

L'intérêt du changement de base est la parcimonie. Pour la rendre concrète, compte combien de déterminants de la base computationnelle sont nécessaires pour capturer 99 % du poids de l'état fondamental, et compare la base des impulsions à la base de position. L'énergie est indépendante de la base, mais pas le nombre de déterminants :

>>> def determinants_for(vec, fraction=0.99):
...     weights = np.sort(np.abs(vec) ** 2)[::-1]
...     return int(np.searchsorted(np.cumsum(weights), fraction) + 1)
...
>>> # the same Hamiltonian in the position basis has an identical spectrum, but a dense ground state
>>> position_linop = ffsim.linear_operator(hamiltonian_position, norb, nelec)
>>> position_energy, position_ground = scipy.sparse.linalg.eigsh(position_linop, k=1, which="SA")
>>>
>>> print(f"position-basis energy: {position_energy[0]:.6f}")
position-basis energy: -13.422492
>>> print(f"determinants for 99% weight: {determinants_for(position_ground[:, 0])} (position)")
determinants for 99% weight: 2563 (position)
>>> print(f"determinants for 99% weight: {determinants_for(ground_state)} (momentum)")
determinants for 99% weight: 21 (momentum)

Les deux opérateurs partagent la même énergie, mais l'état fondamental en base de position est réparti sur des milliers de déterminants, tandis que la base des impulsions le concentre sur environ 20, ce qui explique pourquoi l'échantillonnage SKQD fonctionne dans la base des impulsions et non dans la base de position.

3. Préparer l'état de référence

La référence SKQD $\lvert \psi_0 \rangle$ est la superposition de toutes les excitations des électrons les plus proches du niveau de Fermi vers les modes d'impulsion vides les plus proches. La publication SKQD la construit avec un réseau écrit à la main de XXPlusYYGate. Elle paraît fortement corrélée (elle compte des centaines d'amplitudes non nulles dans la base computationnelle), mais un réseau de rotations à deux modes conservant le nombre de particules est une opération à particule unique (gaussienne fermionique), donc l'état qu'il produit est un unique déterminant de Slater. Tu peux donc le préparer avec une seule porte PrepareSlaterDeterminant, étant donné la rotation à particule unique que le réseau implémente.

Comme chaque XXPlusYYGate agit comme une rotation de Givens sur une paire de modes, cette rotation peut être construite directement dans la représentation à particule unique (une matrice de Givens par porte) au lieu d'émettre à la main les portes à deux qubits. Les portes forment un appareillage en briques de rotations entre plus proches voisins qui s'éventaille depuis le niveau de Fermi en couches de plus en plus larges, répartissant chaque électron proche du niveau de Fermi sur les modes vides situés juste au-dessus :

>>> def fermi_level_rotation(norb, nocc, n_excited=3):
...     """Orbital rotation exciting the near-Fermi electrons into the nearest empty modes."""
...     cos, sin = np.cos(np.pi / 4), np.sin(np.pi / 4)
...     fermi = nocc - 1  # highest occupied mode; the Fermi level sits between it and mode `nocc`
...     # Nearest-neighbor Givens rotations fanning out from the Fermi level in expanding layers:
...     # layer `depth` couples the mode `depth` steps below the level to the one above it, mixing
...     # the `n_excited` highest-occupied modes with the empty modes just above them. A final
...     # rotation carries amplitude up into the topmost empty mode of the window.
...     network = [
...         (lower, lower + 1)
...         for depth in range(n_excited)
...         for lower in range(fermi - depth, fermi + depth + 1, 2)
...     ]
...     network.append((fermi + n_excited, fermi + n_excited + 1))
...     unitary = np.eye(norb, dtype=complex)
...     for i, j in network:
...         givens = np.eye(norb, dtype=complex)
...         givens[i, i] = givens[j, j] = cos
...         givens[i, j], givens[j, i] = sin, -sin
...         unitary = givens @ unitary
...     return unitary
...
>>> reference_rotation = fermi_level_rotation(norb, nelec[0])

Remarque

La rotation ci-dessus est exprimée dans la base des impulsions, en accord avec l'hamiltonien. Les électrons proches du niveau de Fermi sont les modes d'impulsion occupés de plus haute énergie, et le réseau les excite vers les modes vides de plus basse énergie. C'est pourquoi la référence doit être construite après le changement de base de l'étape précédente, et non avant.

La préparation du déterminant de Slater se fait par secteur de spin, puisque chaque secteur a sa propre forme (électrons, orbitales). Une seule porte PrepareSlaterDeterminant (une occupation, les nocc modes les plus bas remplis, avec la rotation à particule unique) prépare un secteur ; la même porte est appliquée une fois sur les modes alpha 0..norb et une fois sur les modes bêta norb..2*norb, puisque les deux secteurs de spin partagent ici la même occupation de référence et la même rotation :

>>> from qiskit_fermions.circuit.library import PrepareSlaterDeterminant
>>>
>>> occupation = [i < nelec[0] for i in range(norb)]
>>> reference_preparation = PrepareSlaterDeterminant(occupation, reference_rotation)

Cette porte déclarative remplace le réseau de XXPlusYYGate codé à la main ; le transpileur génère la synthèse par rotations de Givens.

4. Assembler les circuits de Krylov

Chaque circuit de Krylov est la préparation de l'état de référence suivie d'une évolution sous $e^{-i k \, \Delta t \, H}$. Sépare l'hamiltonien en ses parties à un corps et à deux corps, $H = H_1 + H_2$, et effectue la trotterisation de l'évolution selon cette séparation. C'est le choix décisif pour la profondeur du circuit. L'évolution à un corps $e^{-i t H_1}$ est une opération de fermions libres (gaussienne fermionique), c'est donc une rotation de base à particule unique (une OrbitalRotation d'unitaire $e^{-i t \, \mathbf{h}_1}$) plutôt que quelque chose qu'il faudrait trotteriser terme à terme en Pauli. L'évolution à deux corps $e^{-i t H_2}$ n'agit que sur l'unique terme d'interaction locale. Construire le pas de cette manière permet au transpileur d'émettre un appareillage en briques peu profond de rotations de Givens plus une seule porte d'interaction à deux qubits, la structure peu profonde que la publication SKQD construit à la main.

Nous avons déjà la partie à deux corps. C'est l'opérateur two_body_momentum construit à l'étape 2 (l'unique terme local). L'évolution à un corps $e^{-i t H_1}$ est indépendante du spin, c'est donc une seule OrbitalRotation par secteur, d'unitaire $e^{-i t \, \mathbf{h}_1}$, appliquée une fois à chaque secteur de spin, selon la même structure parallèle que la préparation de l'état de référence.

Chaque circuit de Krylov évolue pendant un temps total $k \, \Delta t$. Cela est réalisé avec un produit de Trotter du second ordre de $k$ pas de taille $\Delta t$, de sorte que l'erreur par pas reste fixe lorsque la dimension de Krylov augmente. Chaque pas place une rotation complète à un corps (une par secteur de spin) entre deux demi-pas de l'évolution à deux corps (peu coûteuse, à terme unique). La dimension de Krylov $D$ est le nombre de ces circuits, avec dim allant de 0 (état de référence seul) à $D - 1$.

>>> import scipy.linalg
>>>
>>> from qiskit_fermions.circuit import FermionicCircuit
>>> from qiskit_fermions.circuit.library import Evolution, OrbitalRotation
>>>
>>> num_modes = 2 * norb
>>> time_step = 0.2
>>>
>>> half_exp_h2 = Evolution(num_modes, two_body_momentum, time_step / 2)
>>> full_exp_h1 = OrbitalRotation(scipy.linalg.expm(-1j * time_step * h1e_momentum))
>>>
>>> def krylov_circuit(dim):
...     # closes over the fermionic gates defined in the outer scope:
...     # `reference_preparation`, `half_exp_h2` and `full_exp_h1`
...     circuit = FermionicCircuit(num_modes)
...     circuit.append(reference_preparation, range(norb))
...     circuit.append(reference_preparation, range(norb, num_modes))
...     for _ in range(dim):  # second-order Trotter product of `dim` steps
...         circuit.append(half_exp_h2, circuit.modes)
...         circuit.append(full_exp_h1, range(norb))
...         circuit.append(full_exp_h1, range(norb, num_modes))
...         circuit.append(half_exp_h2, circuit.modes)
...     return circuit

Remarque

Les circuits construits ici ne comportent pas encore de mesures, car les mêmes circuits fermioniques pilotent aussi la simulation exacte de l'étape 6, et une mesure n'est pas une opération unitaire. L'étape 5 les ajoute avec FermionicCircuit.measure_all(), qui lit le mode i dans le bit classique i.

Préfère mesurer le circuit fermionique plutôt que le circuit transpilé. Une mesure fermionique désigne le mode qu'elle lit, de sorte que le placement et le routage l'accompagnent avec le qubit sur lequel ce mode a fini, et la chaîne de bits échantillonnée est indexée par mode. Mesurer à la place le circuit transpilé te laisse défaire toi-même final_index_layout(), et renvoie silencieusement des résultats permutés si tu l'oublies. Cela repose sur l'encodage de Jordan-Wigner en base d'occupation qu'utilise le préréglage, où chaque mode est porté par un qubit.

5. Transpiler vers des circuits de qubits

La transpilation d'un FermionicCircuit fait correspondre ses portes fermioniques à des portes de qubits. Elle utilise le mapping de Jordan-Wigner jordan_wigner(). Les portes OrbitalRotation se synthétisent en un appareillage en briques de rotations de Givens, et l'Evolution à terme unique se réduit à une seule porte d'interaction à deux qubits.

Comme chaque porte fermionique ici a une synthèse par défaut, tu peux utiliser directement le generate_preset_jw_pass_manager() prêt à l'emploi plutôt que d'assembler les étapes à la main. Les arguments nommés sont transmis à l'étape de qubits ; passe optimization_level=0 pour une image fidèle et non optimisée de la profondeur synthétisée. Générer toute la famille de circuits de Krylov revient alors à une boucle sur la dimension de Krylov. Conserve les FermionicCircuit non transpilés à côté des circuits de qubits transpilés. Les circuits fermioniques pilotent la simulation exacte (vecteur d'état) de l'étape 6, ils doivent donc rester sans mesure ; measure_all(inplace=False) renvoie une copie mesurée à transpiler et laisse l'original intact, et cette copie est ce qu'un backend exécuterait.

>>> from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
>>>
>>> pass_manager = generate_preset_jw_pass_manager(optimization_level=0)
>>>
>>> def krylov_circuits(krylov_dim):
...     fermionic_circuits, circuits = [], []
...     for dim in range(krylov_dim):
...         fermionic_circuit = krylov_circuit(dim)
...         fermionic_circuits.append(fermionic_circuit)
...         measured = fermionic_circuit.measure_all(inplace=False)
...         circuits.append(pass_manager.run(measured))
...     return fermionic_circuits, circuits
...
>>> krylov_dim = 5
>>> fermionic_circuits, circuits = krylov_circuits(krylov_dim)

La dimension de Krylov apparaît directement dans la profondeur à deux qubits des circuits. La préparation de l'état de référence est un coût fixe, et chaque puissance de Krylov supplémentaire ajoute un pas de Trotter du second ordre de l'opérateur d'évolution temporelle.

>>> two_qubit = lambda instr: instr.operation.num_qubits == 2
>>> [circuit.depth(two_qubit) for circuit in circuits]
[4, 12, 21, 30, 39]

Remarque

La séparation de l'hamiltonien maintient ces profondeurs modestes. L'évolution à un corps est une véritable rotation d'orbitales, elle se synthétise donc en un appareillage en briques de XXPlusYYGate sur des qubits adjacents. Chaque demi-pas de l'interaction locale est une unique RZZGate couplant les deux modes de spin de l'impureté ; le terme peu coûteux que le produit du second ordre place de part et d'autre de la rotation. Confier tout l'hamiltonien à une seule porte Evolution trotteriserait au contraire la partie à un corps terme à terme en Pauli, propageant de longues chaînes Z de Jordan-Wigner à travers les couplages impureté-bain (à longue portée) et gonflant la profondeur d'un ordre de grandeur.

Le premier circuit ne prépare que le déterminant de référence, tandis que le dernier comporte le plus de pas de Trotter (quatre dans ce cas) :

>>> circuits[0].draw("mpl", fold=-1, measure_arrows=False)
<Figure size ... with 1 Axes>

>>> circuits[-1].draw("mpl", fold=-1, measure_arrows=False)
<Figure size ... with 1 Axes>

6. Échantillonner des chaînes de bits à partir des circuits de Krylov

Sur du matériel, tu exécuterais les circuits transpilés et tu mesurerais. Voici la version sans bruit (vecteur d'état) de cette étape, qui montre les chaînes de bits que SKQD collecterait. Simule chaque FermionicCircuit non transpilé avec ffsim.apply_unitary() et échantillonne le vecteur d'état obtenu avec ffsim.sample_state_vector().

Une porte PrepareSlaterDeterminant fonctionne en simulation selon le principe valider puis tourner (elle vérifie que l'état entrant occupe son déterminant de référence, puis applique la rotation). Par conséquent, l'entrée est le simple déterminant d'occupation : les nocc modes les plus bas remplis par spin, sans rotation. La porte applique elle-même reference_rotation.

>>> from collections import Counter
>>>
>>> reference = ffsim.slater_determinant(
...     norb, (list(range(nelec[0])), list(range(nelec[1])))
... )
>>>
>>> shots = 1000
>>> counts = Counter()
>>> for dim, fermionic_circuit in enumerate(fermionic_circuits):
...     statevector = ffsim.apply_unitary(reference, fermionic_circuit, norb=norb, nelec=nelec)
...     samples = ffsim.sample_state_vector(
...         statevector, norb=norb, nelec=nelec, shots=shots, seed=dim
...     )
...     counts += Counter(samples)

En regroupant les échantillons des cinq circuits de Krylov, le support est minuscule : quelques centaines de chaînes de bits distinctes sur les $\binom{8}{4}^2 = 4900$ déterminants du secteur (4, 4). Cette concentration, conséquence directe de la parcimonie en base des impulsions établie à l'étape 2, rend tractable la diagonalisation classique qui suit.

Tracer les comptages avec qiskit.visualization.plot_histogram() rend la parcimonie visible. Quelques configurations dominent, en tête desquelles le déterminant de référence.

>>> from qiskit.visualization import plot_histogram
>>>
>>> plot_histogram(dict(counts), number_to_keep=25)
<Figure size ... with 1 Axes>

7. Diagonaliser dans le sous-espace échantillonné

Les chaînes de bits échantillonnées sont l'entrée de la moitié classique de SKQD. L'hamiltonien est projeté sur le sous-espace engendré par les configurations échantillonnées et diagonalisé dans ce sous-espace. Cette étape est confiée au paquet qiskit-addon-sqd, dont diagonalize_fermionic_hamiltonian() exécute la boucle complète de diagonalisation quantique par échantillonnage. Elle construit un sous-espace à partir de lots de configurations échantillonnées, y diagonalise l'hamiltonien, puis affine itérativement le sous-espace grâce à la récupération de configurations : inverser des occupations pour réparer les configurations qui violent la symétrie connue du nombre de particules. La récupération de configurations est conçue pour annuler les dégâts du bruit matériel, qui corrompt les chaînes de bits échantillonnées en les faisant tomber dans un mauvais secteur de nombre de particules ; l'échantillonnage ici est sans bruit, il n'y a donc rien à réparer de ce côté.

Le solveur consomme les comptages échantillonnés sous forme de BitArray, construit directement à partir des counts regroupés de l'étape 6. Passe num_bits explicitement : sans lui, la largeur est déduite de la plus large chaîne de bits échantillonnée, qui est plus étroite que 2 * norb dès qu'aucun échantillon n'occupe le dernier mode, et le solveur rejette alors le tableau. Le solveur diagonalise dans la même base des impulsions que le reste de la construction ; il prend donc directement les intégrales h1e_momentum et h2e_momentum de l'étape 2, sans reconditionnement. L'énergie récupérée coïncide avec la référence exacte à bien moins d'un milli-Hartree près, à partir de quelques centaines de configurations échantillonnées seulement. Comme l'énergie SQD est une borne supérieure variationnelle de l'état fondamental réel, elle se situe toujours au niveau de reference_energy ou au-dessus (toujours disponible depuis l'étape 2) :

>>> from qiskit.primitives import BitArray
>>> from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian
>>>
>>> bit_array = BitArray.from_counts(dict(counts), num_bits=2 * norb)
>>>
>>> # record each iteration's batch results so we can watch the energy converge
>>> result_history = []
>>> result = diagonalize_fermionic_hamiltonian(
...     h1e_momentum,
...     h2e_momentum,
...     bit_array,
...     samples_per_batch=100,
...     norb=norb,
...     nelec=nelec,
...     num_batches=3,
...     max_iterations=3,
...     symmetrize_spin=True,
...     callback=result_history.append,
...     seed=np.random.default_rng(24),
... )
>>> print(f"SQD estimate: {result.energy:.4f} (exact: {reference_energy:.4f})")
SQD estimate: -13.4224 (exact: -13.4225)

Visualise l'exécution : la convergence de l'énergie au fil des itérations, et l'occupation moyenne de chaque orbitale spatiale dans l'état fondamental récupéré. L'énergie par itération est la plus basse parmi les lots de cette itération, et l'occupation additionne les orbital_occupancies sur les deux secteurs de spin. Même sans bruit à corriger, l'énergie diminue d'une itération à l'autre. Chaque itération diagonalise dans num_batches lots de seulement samples_per_batch configurations sous-échantillonnées depuis le pool, et l'affinement itératif dû au report des chaînes de bits améliore régulièrement les configurations qui se retrouvent dans ces lots, de sorte que l'énergie récupérée s'améliore :

>>> import matplotlib.pyplot as plt
>>> from matplotlib.ticker import MaxNLocator
>>>
>>> min_energies = [min(batch, key=lambda r: r.energy).energy for batch in result_history]
>>> occupancies = np.sum(result.orbital_occupancies, axis=0)
>>>
>>> fig, axs = plt.subplots(1, 2, figsize=(12, 5))
>>> _ = axs[0].plot(range(len(min_energies)), min_energies, marker="o", label="SQD energy")
>>> _ = axs[0].axhline(reference_energy, color="gray", linestyle="--", label="exact energy")
>>> _ = axs[0].set_xlabel("iteration")
>>> _ = axs[0].set_ylabel("energy")
>>> _ = axs[0].set_title("Energy vs SQD iteration")
>>> _ = axs[0].xaxis.set_major_locator(MaxNLocator(integer=True))
>>> _ = axs[0].legend()
>>> _ = axs[1].bar(range(norb), occupancies)
>>> _ = axs[1].set_xlabel("spatial orbital")
>>> _ = axs[1].set_ylabel("average occupancy")
>>> _ = axs[1].set_title("Occupancy per spatial orbital")
>>> fig
<Figure size ... with 2 Axes>

Cela clôt la boucle SKQD. Les circuits en base des impulsions des étapes 1 à 5 produisent un échantillon parcimonieux (étape 6), et la diagonalisation par échantillonnage ci-dessus transforme cet échantillon en énergie de l'état fondamental.

Prochaines étapes

Ce guide a exécuté tout le pipeline SKQD sur un simulateur de vecteur d'état sans bruit. Sur du matériel, les circuits transpilés seraient exécutés et mesurés à la place de l'échantillonnage du vecteur d'état de l'étape 6, et les comptages obtenus alimenteraient la diagonalisation de l'étape 7 sans changement. Consulte la documentation Qiskit pour l'exécution de circuits, et les tutoriels de l'addon SQD pour des informations sur le post-traitement par diagonalisation de sous-espace, y compris son comportement sur des échantillons bruités, là où la récupération de configurations fait le plus gros du travail.

Lis le guide sur la relation avec ffsim pour comprendre pourquoi ffsim.apply_unitary() et ffsim.sample_state_vector() fonctionnent nativement sur les circuits fermioniques de ce paquet, et pourquoi le secteur à nombre de particules fixe utilisé à l'étape 6 rend efficace en échantillons la parcimonie en base des impulsions établie à l'étape 2.