Одним из наиболее востребованных приложений квантовых вычислений является моделирование физических процессов, например, моделирование молекул для разработки лекарств и задач многих тел. Эти симуляции обычно включают использование оператора Гамильтона H и времени t, которое используется для эволюции произвольного состояния от |ψ⟩ до e^[iHt]|ψ⟩. Создание схемы, использующей оператор e^[iHt], не является очевидным, но мы можем начать с примера с двумя кубитами.

Моделирование экспоненты σ_i

Предположим, мы хотим смоделировать e^[i σ_z⊗σ_z t]. Схема для этого известна как [1]:

где аргумент вентиля вращения R_z равен −2t.

Показать, что эта схема эквивалентна e^[i σ_z⊗σ_z t], не так уж сложно. Помните, что у нас есть эйлероподобное разложение оператора как e^[i σ_z⊗σ_z t]=cos⁡(t)I+i sin⁡(t) σ_z⊗σ_z (которое работает, поскольку мы имеем, что σ_z^2= Я). Таким образом, мы должны показать, что приведенная выше схема упрощается до того же выражения. Делать это немного долго, но не сложно, как показано ниже:

Это можно легко распространить на любое количество кубитов, т. е. e^[i σ_z⊗σ_z⊗⋯⊗σ_z t] [3]. Теперь, что произойдет, если у нас есть матрица Паули, которая не равна σ_z в показателе степени? Что ж, просто примените самоинверсный вентиль U до и после цепочки CX, такой что σ_k=Uσ_zU, где k∈{x,y}. Эти ворота будут воротами Адамара для σ_x и

для σ_y. И если это личность, вы можете просто оставить этот кубит нетронутым!

Моделирование общих экспонент

Но это просто говорит нам, как моделировать струны Паули, и нам может не быть дан гамильтониан H в этой форме. Однако фактом является то, что для любой матрицы A с размерностью, равной степени двойки (что было бы здесь так, поскольку мы работаем с кубитами), существует некоторое разложение в строку Паули следующего вида [2] :

где

Таким образом, зная, как моделировать струны Паули, мы можем в основном моделировать любой гамильтониан H. Используя струну Паули гамильтониана, мы можем построить разложение Троттера, для которого мы можем создать схему, основанную на простой схеме, которую мы видели выше. Разложение Троттера выглядит следующим образом:

где P_k — тензорное произведение матриц Паули, c_k — коэффициент, соответствующий каждому P_k, t — время, введенное в начале, и N — число, называемое числом Троттера. Чем больше N, тем точнее наше моделирование. Специфика декомпозиции Троттера выходит за рамки этой статьи, но хорошее начало — глава 4.1 [3].

Реализация Qiskit

Имея все необходимые сведения, мы можем приступить к реализации моделирования в Qiskit. Сначала мы собираемся определить простую функцию, которая создает схему цепи, введенную в начале для операторов вида e^[i σ_z⊗σ_z⊗⋯⊗σ_z t].

from qiskit import QuantumCircuit
import numpy as np

def sim_z(t, qc, qubits):
    """
    Add gates to simulate a Pauli Z string exponential
   
    Parameters:
        t: float
            Time parameter to simulate for
        qc: QuantumCircuit
            Circuit to append gates to
        qubits: array
            Array indicating qubits indeces (in order) to 
            append the gates to
    """
    for i in range(len(qubits) - 1):
        qc.cx(qubits[i], qubits[i + 1])
    qc.rz(-2 * t, qubits[-1])
    for i in range(len(qubits) - 1, 0, -1):
        qc.cx(qubits[i - 1], qubits[i])
qc = QuantumCircuit(3)
sim_z(np.pi, qc, [0, 1, 2])
qc.draw()

Что выводит:

q_0: ──■─────────────────────────■──
     ┌─┴─┐                     ┌─┴─┐
q_1: ┤ X ├──■───────────────■──┤ X ├
     └───┘┌─┴─┐┌─────────┐┌─┴─┐└───┘
q_2: ─────┤ X ├┤ Rz(-2π) ├┤ X ├─────
          └───┘└─────────┘└───┘

Теперь мы можем расширить это, чтобы мы могли моделировать любую струну Паули.

def sim_pauli(arr, t, qc, qubits):
    """
    Append gates to simulate any Pauli string
    
    Parameters:
        arr: array
            Array encoding the Pauli string
        t: float
            Time parameter to simulate for
        qc: QuantumCircuit
            Circuit to append gates to
        qubits: array
            Array indicating qubits indeces (in order) to 
            append the gates to
    """
    new_arr = []
    new_qub = []
    for idx in range(len(arr)):
        if arr[idx] != 'I':
            new_arr.append(arr[idx])
            new_qub.append(qubits[idx])

    h_y = 1 / np.sqrt(2) * np.array([[1, -1j], [1j, -1]])
    for i in range(len(new_arr)):
        if new_arr[i] == 'X':
            qc.h(new_qub[i])
        elif new_arr[i] == 'Y':
            qc.unitary(h_y, [new_qub[i]], r'$H_y$')
            
    sim_z(t, qc, new_qub)
    
    for i in range(len(new_arr)):
        if new_arr[i] == 'X':
            qc.h(new_qub[i])
        elif new_arr[i] == 'Y':
            qc.unitary(h_y, [new_qub[i]], r'$H_y$')

И, наконец, мы напишем небольшую функцию, которая принимает гамильтониан H как его разложение Троттера и моделирует его в заданной квантовой схеме.

def sim_ham(hamiltonian, t, qc, qubits, trotter=1):
    """
    Simulates Hamiltonian given as Pauli string
    
    Parameters:
        hamiltonian: dict
            Dictionary encoding the hamiltonian with each
            Pauli product as a key with the coefficient as value
        t: float
            Time parameter to simulate for
        qc: QuantumCircuit
            Circuit to append gates to
        qubits: array
            Array indicating qubits indeces (in order) to 
            append the gates to
        
    """
    temp    = QuantumCircuit(len(qubits))
    delta_t = t / trotter
    for pauli in hamiltonian:
        sim_pauli(pauli, hamiltonian[pauli] * delta_t, temp, range(len(qubits)))
        
    for i in range(trotter):
        qc.compose(temp, qubits, inplace=True)

Анализируя эту функцию, мы видим, что она в основном напоминает эту формулу, которую мы ввели в фоновом режиме:

Во-первых, мы видим, что delta_t соответствует t/N. Затем для каждого члена Паули в строке Паули, представляющей гамильтониан (первая петля соответствует произведению ∏_k, поскольку последовательные элементы соответствуют умножению матриц), мы вызываем нашу функцию sim_pauli, которая добавляет соответствующие элементы, используя time = hamiltonian[pauli] * delta_t, что соответствует c_k т/Н в формуле. Наконец, мы просто добавляем полученную схему N раз во второй цикл, что напоминает возведение в степень.

По сравнению с точным расчетом

Чтобы подвести итог, мы можем увидеть, как наша схема моделирования сравнивается с прямым вычислением e^[iHt]|ψ⟩. Здесь мы проверяем это с помощью гамильтониана H=2XZY+5ZXX+2YXZ, времени 1/(2∗π) и N=50, что дает нам высокую точность.

from qiskit import Aer, execute
from sympy import Matrix

qc = QuantumCircuit(3)
sim_ham({"XZY": 2, "ZXX": 5, "YXZ": 2}, 1 / (2 * np.pi), qc, [0, 1, 2], trotter=50)
qc = qc.reverse_bits()

backend = Aer.get_backend('statevector_simulator')
result  = execute(qc, backend).result()
vec     = result.get_statevector()
vec     = vec / vec[0] # global phase
qvec    = vec / np.linalg.norm(vec) # normalize
Matrix(np.round(qvec, 5))

Теперь, когда у нас есть результат с квантовой схемой, мы можем напрямую вычислить ожидаемый результат и сравнить его.

import scipy

# Start with |0> state
start    = np.zeros(2 ** 3)
start[0] = 1

# Get the matrix corresponding to some Pauli product
# This function can be optimized, but it is left as it
# is for clarity purposes
def get_pauli(string):
    init   = string[0]
    string = string[1:]
    if init == 'X':
        out = np.array([[0, 1], [1, 0]])
    elif init == 'Z':
        out = np.array([[1, 0], [0, -1]])
    elif init == 'Y':
        out = np.array([[0, -1j], [1j, 0]])
    else:
        out = np.eye(2)
    
    for p in string:
        if p == 'X':
            out = np.kron(out, np.array([[0, 1], [1, 0]]))
        elif p == 'Z':
            out = np.kron(out, np.array([[1, 0], [0, -1]]))
        elif p == 'Y':
            out = np.kron(out, np.array([[0, -1j], [1j, 0]]))
        else:
            out = np.kron(out, np.eye(2))

    return out
# Hamiltonian is calculated from the Pauli decomposition
decomp = {"XZY": 2, "ZXX": 5, "YXZ": 2}
H      = np.zeros((2 ** 3, 2 ** 3))
H      = H.astype('complex128')
for pauli in decomp:
    H += decomp[pauli] * get_pauli(pauli)
    
# Hamiltonian is exponentiated and we multiply the starting
# vector by it to get our result
simul  = scipy.linalg.expm(1j * H * (1 / (2 * np.pi)))
vec    = simul @ start
vec    = vec / vec[0] # global phase
cvec   = vec / np.linalg.norm(vec) # normalize
Matrix(np.round(cvec, 5))

Просто сравнив два вектора, мы можем сказать, что наше моделирование было довольно точным! В качестве более строгого математически сравнения мы можем использовать точность состояния, которая также показывает нам, что наша симуляция довольно хороша!

from qiskit.quantum_info import state_fidelity

state_fidelity(qvec, cvec)

Который дает:

0.9999834345555809

Версия этой статьи для ноутбука Jupyter доступна здесь. Если вам понравился этот пост, не забудьте подписаться на Quantum Untangled, чтобы получить больше похожего контента!

Ссылки

  • [1] Подход к гамильтоновому моделированию от Давида Хачатряна на Quantum Computing StackExchange.
  • [2] Можно ли разложить произвольные матрицы по базису Паули? от Danylo Y на Quantum Computing StackExchange.
  • [3] Моделирование гамильтонианов электронной структуры с использованием квантовых компьютеров Джеймсом Д. Уитфилдом, Джейкобом Биамонте и Аланом Аспуру-Гузиком.