Sintetizar uma porta Evolution no espaço bosônico

A evolução temporal $e^{-itH}$ de um Hamiltoniano que é uma soma de partes não comutativas não pode ser implementada exatamente porta a porta. Uma fórmula de produto a aproxima por uma sequência de evoluções curtas sob cada uma das partes, e o erro da aproximação é chamado de erro de Trotter.

Você pode aplicar uma fórmula de produto em dois lugares, e eles não são intercambiáveis:

1. No espaço bosônico, antes de qualquer truncamento ou codificação ser escolhido. A fórmula enxerga os grupos do operador (veja Agrupar termos de operadores: use a estrutura do operador) e pode intercalá-los.

2. No espaço de qubits, dentro da PauliEvolutionGate do Qiskit, depois que o operador foi mapeado. A fórmula enxerga um operador de Pauli por grupo e só pode subdividir os termos dentro dele.

Somente a primeira pode reduzir o erro entre grupos que não comutam, que costuma ser o dominante. Este guia mostra como usá-la.

Importante

Os dois passos são opcionais. Uma porta Evolution que nunca é decomposta no espaço bosônico chega inteira à etapa de bóson para qubit e é mapeada e sintetizada ali. Isso não gera nenhum erro no nível bosônico, apenas o que a fórmula no nível de qubits introduzir.

Configuração

Uma cadeia de Bose-Hubbard de dois sítios, o menor exemplo em que os grupos não comutam (veja Simular a dinâmica de Bose-Hubbard com a escolha de uma codificação para o modelo completo). Ela é agrupada em um grupo hermitiano para a ligação e um grupo diagonal para a interação no sítio:

>>> 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 um erro, você precisa de uma referência. A evolução exata do Hamiltoniano como mapeado no mesmo truncamento é a adequada, porque o truncamento é uma aproximação separada:

>>> 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 o circuito convertido sem decompor suas instâncias de PauliEvolutionGate. Statevector aplica cada uma como a exponencial exata do seu grupo, então a única coisa que sobra é o produto de exponenciais de grupos que a síntese bosônica escolheu.

O padrão

Sem o argumento synthesis, uma porta Evolution é decomposta com BosonicLieTrotter, a fórmula de primeira ordem: um fator por grupo, cada um pelo tempo inteiro.

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

Sem decompor, a porta inteira é entregue à etapa de qubits e não há nenhum erro de Trotter bosônico:

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

Decomposta, aparece um erro genuíno, porque os dois grupos não comutam:

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

Reduza o erro

Passe um método de síntese diferente para a porta. Uma ordem maior compõe os fatores de forma simétrica, de modo que os termos de erro dominantes se cancelam; mais repetições dividem o tempo em passos menores.

>>> 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 de duas unidades na ordem compõe cinco varreduras da ordem abaixo (o crescimento $5^{(k-2)/2}$ da recursão de Suzuki, aplicado aos $2n - 1$ fatores da ordem 2), que é a contagem na última linha, marcada com uma estrela. As varreduras se encontram no mesmo fator, então as portas vizinhas evoluem sob um mesmo operador e são fundidas em uma única porta com o tempo somado. Isso economiza uma porta por junção e é exato. Passe merge_adjacent=False para manter a sequência de livro-texto. De qualquer forma, uma ordem maior compra um erro 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

As repetições trocam profundidade por precisão a uma taxa diferente. O erro de primeira ordem cai linearmente com o tamanho do passo, e o de segunda ordem, 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

Qual é a melhor troca depende do operador. Observe que um operador com um único grupo é exato em qualquer ordem, porque um fator isolado comuta consigo mesmo, e a fórmula então se reduz a uma porta em vez de dividi-la novamente.

Meça o custo

Um erro menor custa mais fatores, e cada fator custa portas emaranhadoras depois que a etapa de qubits o sintetiza. Quantas depende da codificação, então conte-as pelo pipeline do preset, com a base fixa e sem otimização no nível de qubits (optimization_level=0), para ver o preço 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

As contagens exatas dependem de como o Qiskit sintetiza as portas de evolução, o que muda entre versões, então o gráfico as mostra em vez de uma tabela:

Cada grupo de barras é um método: Lie-Trotter, Suzuki de ordem 2 e Suzuki de ordem 4. A ordem compra precisão a um preço que cresce pelo fator da recursão, e a codificação unária paga mais que o código de Gray em todos os métodos, porque cada um de seus fatores é mapeado para mais portas emaranhadoras. Escolha o método pela precisão de que você precisa e a codificação pelo hardware, e mantenha as duas escolhas separadas: o circuito bosônico acima nunca muda.

Aplique um método a um circuito inteiro

Para definir synthesis em cada porta, você precisa propagar a escolha por tudo que constrói uma. O passo B2BSynthesis faz a escolha uma única vez para um pipeline. Seu config nomeia o método para cada tipo de porta, e ele tanto seleciona o método quanto expande as portas:

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

Um método também pode ser indicado pelo nome, com os argumentos de seu __init__, em vez de ser construído: B2BSynthesis({"Evolution": ("SuzukiTrotter", (), {"order": 2})}). Os nomes disponíveis são os do B2BSynthesisPluginManager, e um método que não está registrado em lugar nenhum pode ser entregue como uma instância, como acima.

Com apply=False, o passo apenas seleciona o método e deixa a expansão para o ponto em que ela for desejada, por exemplo, para inspecionar o circuito antes. Uma porta atômica nunca é tocada, nem sequer marcada:

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

Uma porta não expandida chega inteira à etapa de bóson para qubit, onde é mapeada sem que seu synthesis seja consultado, então todos os métodos produzem então uma saída idêntica.

Um filter opcional o restringe às portas que se beneficiam. Evoluções cujos grupos todos comutam, como um Hamiltoniano feito apenas de operadores de número, são sintetizadas exatamente em qualquer ordem, então uma ordem maior apenas acrescentaria profundidade ali.

Restrinja o passo às portas que se beneficiam

Um operador puramente diagonal torna o ponto concreto, porque operadores de número comutam entre si. Seus grupos são exatos em primeira ordem, e a segunda ordem gasta portas emaranhadoras à toa:

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

Os dois circuitos são exatos, então o segundo gasta mais portas emaranhadoras sem nenhum ganho. (A fórmula simétrica de ordem dois aplica o primeiro grupo duas vezes, com metade do tempo cada, e por isso o custo depende de qual grupo vem primeiro.) O filter recebe um predicado sobre o DAGOpNode e deixa intocadas as portas rejeitadas, de modo que um único passo pode elevar a ordem das portas que se beneficiam e manter o restante na primeira ordem. O predicado abaixo aceita uma porta somente quando seu operador move bósons entre modos, que é quando os grupos deixam de comutar:

>>> 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 você executar o passo novamente, nada muda: os fatores que ele emite são atômicos.

Fatores atômicos

Uma Evolution marcada como atomic não tem definição, então decompose() a deixa no lugar por mais vezes que seja aplicado. Uma síntese marca dessa forma os fatores que emite, e o BosonicQDrift também faz isso para os termos que amostra.

Isso é mais do que organização. Dividir um grupo hermitiano ainda mais, termo a termo, produz fatores que individualmente não são hermitianos, e a exponencial de um operador não hermitiano não é unitária:

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

Cuidado

Cada fator precisa ser hermitiano por si só. Um operador hermitiano dividido em grupos pode produzir grupos não hermitianos, como um termo de salto separado de seu adjunto. Verifique com groups_are_hermitian(). Um operador sem grupos é dividido termo a termo, e a etapa de qubits então rejeita os fatores com Operator contains complex coefficients. A simetrização de uma ordem par cancela parcialmente o erro de um fator não hermitiano, então um operador agrupado incorretamente pode parecer melhor na ordem 2 do que na ordem 1 e ainda assim estar errado.

Repita um passo curto

Como alternativa no nível do circuito, construa um passo de Trotter curto e repita-o. repeat() achata as repetições em um único circuito e deixa as fronteiras adjacentes, de modo que os passos possam atuar através delas:

>>> 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 para manter as repetições separadas. Uma barreira é mapeada para cada qubit dos registradores que ela toca, porque uma codificação precisa de vários qubits por modo e um modo não tem uma imagem por qubit própria:

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

Nota

repeat não tem relação com o reps de um método de síntese. Uma síntese subdivide uma única evolução em mais fatores, e mais curtos. Repetir um circuito que contém uma evolução apenas multiplica seu tempo total.

Combine com uma fórmula no nível de qubits

Os dois níveis se compõem e, onde ambos aproximam, a precisão do resultado é governada pelo mais fraco dos dois. Elevar a ordem bosônica enquanto a fórmula do lado dos qubits continua de primeira ordem traz pouco ganho. Veja MapperFnEvolutionSynthesis para a fórmula do lado dos qubits. Ela mantém a ordem dos termos de Pauli que o mapeador escolheu, e assim group_wise() pode emitir os termos de cada grupo adjacentes entre si, e simplify() pode fixar uma ordem canônica.

Próximos passos

• Transpilar circuitos bosônicos explica as etapas que se seguem à síntese bosônica e onde B2BSynthesis se encaixa entre elas.

• Agrupar termos de operadores: use a estrutura do operador explica como atribuir e verificar os grupos que uma fórmula de produto intercala.

• Simular a dinâmica de Bose-Hubbard com a escolha de uma codificação aplica o esquema de repetição a um exemplo completo.