Ruotare i modi bosonici

Una ModeRotation mescola i modi bosonici con una matrice unitaria $U$, come fa una rete di beamsplitter e sfasatori. Viene detta passiva perché conserva il numero totale di bosoni: sposta i bosoni tra i modi e non ne crea né ne distrugge mai nessuno. (Lo squeezing e il displacement sono le operazioni attive, e una rotazione di modi non può esprimerle.) Mappa gli operatori di creazione come

$$b^\dagger_i \mapsto \sum_j U_{ji} \, b^\dagger_j \, ,$$

che sull'intero spazio di Fock è l'unitario a molti corpi

$$\exp\left(\sum_{ij} \log(U)_{ij} \, b^\dagger_i b_j\right) \, .$$

Questo è il gate alla base di un interferometro lineare-ottico: invii stati di Fock, ruoti i modi e leggi come i bosoni si ridistribuiscono. Come ogni gate di un BosonicCircuit (vedi Lavorare con i circuiti bosonici), non porta con sé né troncamento né codifica. Entrambi intervengono al momento del transpile, e questa guida mostra perché il troncamento richiede particolare attenzione per una rotazione.

Costruire una rotazione

Il costruttore accetta la matrice unitaria $M \times M$. Il gate agisce su $M$ modi. La classe non verifica che la matrice sia unitaria, quindi è una tua responsabilità.

>>> import numpy as np
>>> from qiskit_bosons.circuit.library import ModeRotation
>>>
>>> unitary = np.array([[1, -1], [1, 1]]) / np.sqrt(2)
>>> rotation = ModeRotation(unitary)
>>> rotation.num_modes  # the gate acts on two modes
2

Le rotazioni più comuni hanno costruttori dedicati. ModeRotation.beamsplitter() prende un angolo di mescolamento $\theta$ e una fase opzionale $\phi$, e $\theta = \pi/4$ dà il beamsplitter bilanciato 50:50. ModeRotation.phase_shift() prende una fase per modo.

>>> beamsplitter = ModeRotation.beamsplitter(np.pi / 4)
>>> np.round(beamsplitter.rotation_unitary.real, 4)
array([[ 0.7071, -0.7071],
       [ 0.7071,  0.7071]])
>>> bool(np.allclose(beamsplitter.rotation_unitary, unitary))
True
>>> shift = ModeRotation.phase_shift([0.0, np.pi / 2])
>>> np.round(np.diag(shift.rotation_unitary), 4)
array([1.+0.j, 0.+1.j])

Puoi anche partire dal generatore anti-hermitiano $K$ e lasciare che il gate lo esponenzi, il che dà la rotazione $\exp(K)$. Il generatore si recupera con generator().

>>> generator = np.array([[0, -0.3], [0.3, 0]], dtype=complex)
>>> small = ModeRotation.from_generator(generator)
>>> np.round(small.rotation_unitary.real, 4)
array([[ 0.9553, -0.2955],
       [ 0.2955,  0.9553]])
>>> bool(np.allclose(small.generator(), generator))
True

Il ponte verso l'algebra degli operatori è to_operator(), che restituisce il generatore hermitiano $H = \sum_{ij} \mathrm{i} \log(U)_{ij} \, b^\dagger_i b_j$ come BosonOperator, in modo che il gate sia $e^{-\mathrm{i} H}$. Per il beamsplitter 50:50, si tratta di una coppia di termini di hopping con coefficienti immaginari opposti. L'operatore è hermitiano, come richiede un gate Evolution, e conserva il numero di particelle.

>>> operator = beamsplitter.to_operator()
>>> # The stored order of the terms is not specified, so sort them before printing.
>>> sorted((term[0], round(term[1].imag, 4)) for term in operator.iter_terms())
[([(True, 0), (False, 1)], -0.7854), ([(True, 1), (False, 0)], 0.7854)]
>>> operator.is_hermitian()
True
>>> operator.conserves_particle_number()
True

Ogni termine si legge (is_creation, mode) da sinistra a destra, quindi il secondo è $\mathrm{i} \tfrac{\pi}{4} \, b^\dagger_1 b_0$.

Agire sugli stati di Fock

Una rotazione sposta i bosoni tra i modi, quindi il suo effetto si vede meglio sugli stati di Fock. L'effetto Hong-Ou-Mandel è l'esempio standard. Invia un bosone in ciascun modo di ingresso di un beamsplitter bilanciato. I due cammini che lasciano i bosoni in modi di uscita diversi interferiscono in modo distruttivo, quindi lo stato $\lvert 1, 1 \rangle$ scompare e i bosoni escono sempre insieme, in $\lvert 2, 0 \rangle$ o $\lvert 0, 2 \rangle$, ciascuno con probabilità un mezzo.

Costruisci il circuito sui modi bosonici e abbassalo con il pass manager preimpostato. Il LevelSpec fissa a questo punto il troncamento e la codifica. Un cutoff di max_occupation=2 mantiene i livelli 0, 1 e 2 di ogni modo.

>>> from qiskit_bosons.circuit import BosonicCircuit
>>> from qiskit_bosons.circuit.library import InitializeModes
>>> from qiskit_bosons.mappers import LevelEncoding, LevelSpec
>>> from qiskit_bosons.transpiler.presets import generate_preset_occupation_pass_manager
>>>
>>> circuit = BosonicCircuit(2)
>>> circuit.append(InitializeModes([1, 1]), circuit.modes)
>>> circuit.append(ModeRotation.beamsplitter(np.pi / 4), circuit.modes)

La pipeline restituisce un normale circuito di qubit. Leggi le probabilità degli esiti dallo stato finale, qui con lo Statevector di Qiskit, che a questa dimensione è esatto. Funziona qualsiasi simulatore che restituisca ampiezze o conteggi. Uno stato di Fock $\lvert n_0, n_1 \rangle$ è un singolo stato di base computazionale, il cui indice è dato dagli schemi di bit di LevelEncoding.level_to_bits(); i qubit del secondo modo seguono quelli del primo.

>>> from qiskit.quantum_info import Statevector
>>> from qiskit.transpiler import PassManager
>>>
>>> def fock_probabilities(circuit, max_occupation):
...     """Maps each pair of occupations to the probability of that outcome."""
...     encoding = LevelEncoding.UNARY_UNPROJECTED
...     spec = LevelSpec(max_occupation, encoding)
...     pass_manager = generate_preset_occupation_pass_manager(spec)
...     pass_manager.qubit = PassManager()  # keep the evolutions as the pipeline emitted them
...     qubit_circuit = pass_manager.run(circuit)
...     state = Statevector(qubit_circuit)
...     width = encoding.num_qubits(max_occupation + 1)
...     probabilities = {}
...     for n0 in range(max_occupation + 1):
...         for n1 in range(max_occupation + 1):
...             index = encoding.level_to_bits(n0) | (encoding.level_to_bits(n1) << width)
...             probability = round(float(abs(state.data[index]) ** 2), 4)
...             if probability > 0:
...                 probabilities[(n0, n1)] = probability
...     return probabilities
>>>
>>> fock_probabilities(circuit, max_occupation=2)
{(0, 2): 0.5, (2, 0): 0.5}

L'esito $\lvert 1, 1 \rangle$ ha probabilità zero, e i due esiti raggruppati si dividono equamente la probabilità.

Scegliere il cutoff

Il risultato precedente richiede max_occupation=2. Esegui lo stesso circuito con un cutoff di 1. Basta a contenere lo stato di ingresso $\lvert 1, 1 \rangle$, ma non lo stato che la rotazione ne produce, che ha peso su $\lvert 2, 0 \rangle$ e $\lvert 0, 2 \rangle$:

>>> fock_probabilities(circuit, max_occupation=1)
{(1, 1): 1.0}

Il circuito viene eseguito senza alcun errore, eppure la risposta è sbagliata: lo stato resta in $\lvert 1, 1 \rangle$ con probabilità uno. Una rotazione mescola settori di occupazione all'interno di un modo. Due bosoni possono finire nello stesso modo di uscita, quindi lo stato esatto ha peso sull'occupazione 2 in un modo che era partito con un bosone. Con un cutoff di 1 quel livello non esiste, e il generatore troncato non può spostare lo stato lì.

Un limite sull'occupazione dello stato iniziale non è quindi sufficiente. Il cutoff deve coprire il numero totale di bosoni che può raggiungere un qualsiasi singolo modo della rotazione. Per $N$ bosoni in totale, la scelta sicura è un cutoff pari a $N$.

Lo stesso ragionamento vale per ogni modo separatamente, quindi una combinazione di cutoff diversi per modo (vedi Descrivere il troncamento con un LevelSpec) rende facile mappare una rotazione in modo errato. Un modo che può ricevere bosoni richiede un cutoff che copra ciò che può ricevere, non solo ciò che contiene inizialmente. Il controllo LevelSpec.truncation() confronta una specifica con il limite inferiore di occupazione di un operatore. Non può segnalare il caso visto sopra. Ogni termine del generatore del beamsplitter contiene un singolo operatore di creazione, quindi il limite è uno e il controllo non riporta alcuna carenza con un cutoff di 1, sebbene il risultato lì sia sbagliato:

>>> LevelSpec(1, "unary_unprojected").truncation(beamsplitter.to_operator())
[0, 0]

Un'ulteriore precauzione deriva dal modo in cui il transpiler costruisce una rotazione (sezione successiva). La applica come una sequenza di fattori a due modi, e la documentazione del plugin avverte che un singolo fattore può popolare un livello che la rotazione completa non raggiunge mai. La sequenza troncata e la singola evoluzione troncata sono allora approssimazioni diverse dello stesso unitario. Una ricerca su 1.080 rotazioni casuali a tre modi, con un bosone in ciascuno di due modi di ingresso e nessuna doppia occupazione in uscita, non ha rilevato un caso simile, quindi consideralo una precauzione e non un fallimento comune. Fai convergere il cutoff come per qualsiasi altro circuito, aumentandolo finché il risultato non smette di cambiare.

Come il transpiler sintetizza una rotazione

La fase optimization della pipeline preimpostata scompone una ModeRotation prima nello spazio bosonico, con il pass B2BSynthesis e il metodo BosonicClementsDecomposition (registrato come ModeRotation.Clements). Questa è una scelta del preset e non un requisito del framework: una sezione successiva abbassa lo stesso gate senza alcuna decomposizione. La decomposizione non esponenzia il generatore completo in un solo passo. Prima fattorizza l'unitario in al più $M(M-1)/2$ rotazioni tra modi adiacenti più $M$ fasi a singolo modo, seguendo la maglia rettangolare di Clements et al. [1]. Ogni fattore diventa un gate Evolution atomico, che la fase di conversione da bosone a qubit mappa come qualsiasi altro.

La fattorizzazione vera e propria è givens_decomposition(), che la decomposizione richiama. La mappa da un'unitaria alla rotazione a molti corpi è un omomorfismo di gruppi, quindi il prodotto delle rotazioni dei fattori è esattamente uguale alla rotazione completa sullo spazio non troncato. La mesh quindi non aggiunge alcun errore proprio. Il suo vantaggio è la località. Ogni evoluzione emessa coinvolge solo i termini di Pauli dei due modi che mescola (o di un solo modo per una fase), invece dei termini su tutti i modi.

Esegui il transpile di una rotazione casuale a tre modi e ispeziona ciò che la pipeline ha emesso. Una unitaria generica $3 \times 3$ richiede al massimo tre rotazioni, più tre fasi, e ogni fattore è un proprio PauliEvolutionGate. Tutti sono dichiarati sull'intero registro, ma gli operatori al loro interno agiscono solo su pochi qubit:

>>> from qiskit.quantum_info import random_unitary
>>>
>>> random = np.asarray(random_unitary(3, seed=7).data)
>>> spec = LevelSpec(2, "unary_unprojected")
>>> rotation_circuit = BosonicCircuit(3)
>>> rotation_circuit.append(InitializeModes([1, 0, 0]), rotation_circuit.modes)
>>> rotation_circuit.append(ModeRotation(random), rotation_circuit.modes)
>>> pass_manager = generate_preset_occupation_pass_manager(spec)
>>> pass_manager.qubit = PassManager()
>>> qubit_circuit = pass_manager.run(rotation_circuit)
>>> dict(qubit_circuit.count_ops())
{'PauliEvolution': 6, 'x': 3}
>>> for instruction in qubit_circuit.data:
...     if instruction.operation.name == "PauliEvolution":
...         terms = instruction.operation.operator.to_sparse_list()
...         print(len(terms), sorted({qubit for _, indices, _ in terms for qubit in indices}))
64 [0, 1, 2, 3, 4, 5]
64 [3, 4, 5, 6, 7, 8]
64 [0, 1, 2, 3, 4, 5]
2 [1, 2]
2 [4, 5]
2 [7, 8]

Con tre livelli per modo, ogni modo occupa tre qubit. Un beamsplitter tra due modi coinvolge i sei qubit di quei modi, e una fase coinvolge i qubit di un solo modo.

I fattori sono esatti solo quando ogni evoluzione è sintetizzata bene. I termini del generatore di un beamsplitter non commutano, quindi la formula di prodotto predefinita a un solo passo lascia un errore in ogni fattore. Passa una formula di prodotto con più ripetizioni al plugin Evolution che mappa i fattori per ridurlo:

>>> from qiskit.passmanager import MultiStagePassManager
>>> from qiskit.synthesis import LieTrotter
>>> from qiskit_bosons.circuit.library.synthesis import BosonicClementsDecomposition
>>> from qiskit_bosons.mappers.library import occupation_encoding
>>> from qiskit_bosons.transpiler import (
...     BosonicCircuitToDAG,
...     BosonicPassManager,
...     QuantumDAGToCircuit,
... )
>>> from qiskit_bosons.transpiler.passes import (
...     B2BSynthesis,
...     B2QSynthesis,
...     MapperFnEvolutionSynthesis,
...     OccupationB2QLayout,
...     OccupationInitializeModesSynthesis,
... )
>>>
>>> def with_product_formula(product_formula):
...     synthesis = B2QSynthesis()
...     synthesis.methods["Evolution"] = MapperFnEvolutionSynthesis(
...         lambda operator, num_qubits: occupation_encoding(operator, spec, num_qubits),
...         product_formula,
...     )
...     synthesis.methods["InitializeModes"] = OccupationInitializeModesSynthesis(spec)
...     return MultiStagePassManager(
...         input=BosonicCircuitToDAG(),
...         optimization=BosonicPassManager(
...             B2BSynthesis({"ModeRotation": BosonicClementsDecomposition()})
...         ),
...         layout=OccupationB2QLayout(spec),
...         synthesis=synthesis,
...         output=QuantumDAGToCircuit(),
...     )

Un singolo bosone che entra nel modo 0 esce nel modo $j$ con probabilità $|U_{j0}|^2$. Confronta il circuito sintetizzato con questo risultato esatto per numeri crescenti di ripetizioni. Decomponi prima il circuito in modo che la formula di prodotto venga applicata; altrimenti un simulatore esatto come Statevector applica il gate come esponenziale completa e nasconde l'errore.

>>> encoding = LevelEncoding.UNARY_UNPROJECTED
>>> width = encoding.num_qubits(3)
>>> exact = np.abs(random[:, 0]) ** 2
>>>
>>> def one_boson_probabilities(qubit_circuit):
...     """Probability that the single boson sits in mode j, for each j."""
...     state = Statevector(qubit_circuit)
...     probabilities = []
...     for j in range(3):
...         occupations = [1 if mode == j else 0 for mode in range(3)]
...         index = sum(
...             encoding.level_to_bits(n) << (width * mode) for mode, n in enumerate(occupations)
...         )
...         probabilities.append(abs(state.data[index]) ** 2)
...     return np.array(probabilities)
>>>
>>> for repetitions in (1, 4, 16):
...     lowered = with_product_formula(LieTrotter(reps=repetitions)).run(rotation_circuit)
...     error = np.max(np.abs(one_boson_probabilities(lowered.decompose(reps=2)) - exact))
...     print(f"{repetitions:2d} repetitions: error {error:.1e}")
 1 repetitions: error 1.2e-03
 4 repetitions: error 7.6e-05
16 repetitions: error 4.7e-06

L'errore diminuisce di un fattore circa quattro per ogni quadruplicazione delle ripetizioni, come atteso per una formula del primo ordine. Il costo cresce in proporzione, perché ogni fattore viene ripetuto.

Abbassare una rotazione senza decomporla

Una decomposizione nello spazio bosonico è un modo per abbassare una rotazione. (La libreria non include un altro modo, ma la fase da bosone a qubit non dipende da esso.) B2QSynthesis cerca un plugin in base al nome del gate, e un ModeRotation lasciato intero lo raggiunge come qualsiasi altro gate. Un plugin per esso mappa in un colpo solo il generatore di ModeRotation.to_operator() ed emette un solo PauliEvolutionGate a tempo unitario. Tutto qui:

>>> from qiskit.circuit.library import PauliEvolutionGate
>>> from qiskit_bosons.transpiler.passes.utils import map_node_single_register
>>>
>>> class DirectModeRotationSynthesis:
...     """Maps the generator of the whole rotation, with no bosonic decomposition."""
...
...     def __init__(self, mapper_fn):
...         self.mapper_fn = mapper_fn
...
...     def run(self, in_node, out_dag, *, b2q_layout):
...         modes, qreg = map_node_single_register(in_node, b2q_layout)
...         generator = in_node.op.to_operator().relabel_modes(list(modes))
...         observable = self.mapper_fn(generator, len(qreg))
...         out_dag.apply_operation_back(PauliEvolutionGate(observable, time=1.0), tuple(qreg))

Registralo per il nome del gate in un pass manager che non lo decomponga:

>>> def mapper_fn(operator, num_qubits):
...     return occupation_encoding(operator, spec, num_qubits)
>>>
>>> synthesis = B2QSynthesis()
>>> synthesis.methods["Evolution"] = MapperFnEvolutionSynthesis(mapper_fn)
>>> synthesis.methods["InitializeModes"] = OccupationInitializeModesSynthesis(spec)
>>> synthesis.methods["ModeRotation"] = DirectModeRotationSynthesis(mapper_fn)
>>> direct = MultiStagePassManager(
...     input=BosonicCircuitToDAG(),
...     layout=OccupationB2QLayout(spec),
...     synthesis=synthesis,
...     output=QuantumDAGToCircuit(),
... )
>>> direct_circuit = direct.run(rotation_circuit)
>>> dict(direct_circuit.count_ops())
{'x': 3, 'PauliEvolution': 1}
>>> terms = direct_circuit.data[-1].operation.operator.to_sparse_list()
>>> len(terms), sorted({qubit for _, indices, _ in terms for qubit in indices})
(198, [0, 1, 2, 3, 4, 5, 6, 7, 8])

La rotazione è ora un'unica evoluzione su tutti e nove i qubit, mentre la mesh ne ha sei, nessuna delle quali si estende su più di sei. Le due non sono la stessa approssimazione. Il troncamento viene applicato una sola volta al generatore dell'intera rotazione qui, e a ogni fattore separatamente nella mesh (da cui l'avvertenza sul cutoff sopra). Il percorso diretto è uguale a occupation_encoding() dell'operatore a qualsiasi cutoff, quindi è il riferimento con cui confrontare una mesh. Il suo costo è la località: i suoi termini si estendono su tutti i modi, e una formula di prodotto ne ha di più da suddividere.

Nessuna delle due è migliore in generale, ed è per questo che conviene lasciare aperta la scelta. Scrivi un plugin come questo per confrontare gli abbassamenti, per usare una rotazione specifica della codifica, o per passare il gate a qualcosa che lo comprenda nativamente.

Passi successivi

• Esegui il transpile di circuiti bosonici spiega le fasi eseguite dal pass manager predefinito e come sostituire il plugin che gestisce un gate.

• Mappa le rappresentazioni degli operatori tratta il troncamento, le codifiche e il LevelSpec che fissa il layout di ogni modo.

• Lavora con i circuiti bosonici mostra come assemblare e combinare circuiti bosonici.

• La documentazione API elenca qiskit_bosons.circuit.library per i gate, qiskit_bosons.linalg per la decomposizione e qiskit_bosons.transpiler.passes per i plugin di sintesi.

[1]

W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S. Kolthammer, and I. A. Walmsley, Optimal design for universal multiport interferometers, Optica 3, 1460 (2016), arXiv:1603.08788.