Екземпляри та розширення

У цьому розділі розглянемо кілька квантових варіаційних алгоритмів, зокрема:

• Variational Quantum Eigensolver (VQE)
• Subspace Search VQE (SSVQE)
• Variational Quantum Deflation (VQD)
• Quantum Sampling Regression (QSR)

Вивчаючи ці алгоритми, ти ознайомишся з кількома ідеями проєктування, які можна включати до власних варіаційних алгоритмів: ваги, штрафи, надвибірка та недовибірка. Заохочуємо тебе експериментувати з цими концепціями та ділитися результатами з спільнотою.

Фреймворк шаблонів Qiskit застосовується до всіх цих алгоритмів — але кроки ми явно позначатимемо лише в першому прикладі.

Variational Quantum Eigensolver (VQE)

VQE — один із найбільш широко використовуваних варіаційних квантових алгоритмів, що задає шаблон, на основі якого будуються інші алгоритми.

Крок 1: Відображення класичних вхідних даних на квантову задачу

Теоретична схема

Схема VQE проста:

• Підготуй референсні оператори $U_R$
  • Починаємо зі стану $|0\rangle$ і переходимо до референсного стану $|\rho\rangle$
• Застосуй варіаційну форму $U_V(\vec\theta_{i,j})$, щоб створити анзац $U_A(\vec\theta_{i,j})$
  • Переходимо зі стану $|\rho\rangle$ до $U_V(\vec\theta_{i,j})|\rho\rangle = |\psi(\vec\theta_{i,j})\rangle$
• Виконай ініціалізацію при $i=0$, якщо маєш подібну задачу (зазвичай знайдену через класичне моделювання або вибірку)
  • Кожен оптимізатор ініціалізуватиметься по-різному, породжуючи початковий набір векторів параметрів $\Theta_0 := \\{ {\vec\theta_{0,j} | j \in \mathcal{J}_\text{opt}^0} \\}$ (наприклад, з початкової точки $\vec\theta_0$).
• Обчисли функцію вартості $C(\vec\theta_{i,j}) := \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ для всіх підготовлених станів на квантовому комп'ютері.
• Використай класичний оптимізатор для вибору наступного набору параметрів $\Theta_{i+1}$.
• Повторюй процес до досягнення збіжності.

Це простий цикл класичної оптимізації, де ми обчислюємо функцію вартості. Деякі оптимізатори можуть вимагати кількох обчислень для розрахунку градієнта, визначення наступної ітерації або оцінки збіжності.

Ось приклад для такого спостережуваного оператора:

$$\hat{O}_1 = 2 II - 2 XX + 3 YY - 3 ZZ,$$

Реалізація

# Added by doQumentation — required packages for this notebook
!pip install -q numpy qiskit qiskit-ibm-runtime scipy

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit.circuit.library import TwoLocal
import numpy as np

theta_list = (2 * np.pi * np.random.rand(1, 8)).tolist()
observable = SparsePauliOp.from_list([("II", 2), ("XX", -2), ("YY", 3), ("ZZ", -3)])

reference_circuit = QuantumCircuit(2)
reference_circuit.x(0)

variational_form = TwoLocal(
    2,
    rotation_blocks=["rz", "ry"],
    entanglement_blocks="cx",
    entanglement="linear",
    reps=1,
)

ansatz = reference_circuit.compose(variational_form)

ansatz.decompose().draw("mpl")

def cost_func_vqe(parameters, ansatz, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """

    estimator_job = estimator.run([(ansatz, hamiltonian, [parameters])])
    estimator_result = estimator_job.result()[0]

    cost = estimator_result.data.evs[0]
    return cost

from qiskit.primitives import StatevectorEstimator

estimator = StatevectorEstimator()

Ми можемо використати цю функцію вартості для обчислення оптимальних параметрів:

# SciPy minimizer routine
from scipy.optimize import minimize

x0 = np.ones(8)

result = minimize(
    cost_func_vqe, x0, args=(ansatz, observable, estimator), method="COBYLA"
)

result

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -5.999999982445723
       x: [ 1.741e+00  9.606e-01  1.571e+00  2.115e-05  1.899e+00
            1.243e+00  6.063e-01  6.063e-01]
    nfev: 136
   maxcv: 0.0

Крок 2: Оптимізація задачі для квантового виконання

Обираємо найменш завантажений Backend і імпортуємо необхідні компоненти з qiskit_ibm_runtime.

from qiskit_ibm_runtime import SamplerV2 as Sampler
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit_ibm_runtime import Session, EstimatorOptions
from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
print(backend)

<IBMBackend('ibm_brisbane')>

Транспілюємо Circuit за допомогою менеджера попередніх проходів з рівнем оптимізації 3 і застосовуємо відповідне розміщення до спостережуваного оператора.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
isa_ansatz = pm.run(ansatz)
isa_observable = observable.apply_layout(layout=isa_ansatz.layout)

Крок 3: Виконання за допомогою примітивів Qiskit Runtime

Тепер ми готові запустити обчислення на залізі IBM Quantum®. Оскільки мінімізація функції вартості є дуже ітеративною, ми відкриємо сесію Runtime. Завдяки цьому нам потрібно буде чекати в черзі лише один раз. Після того як завдання почне виконуватися, кожна ітерація з оновленими параметрами запускатиметься негайно.

x0 = np.ones(8)

estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    result = minimize(
        cost_func_vqe,
        x0,
        args=(isa_ansatz, isa_observable, estimator),
        method="COBYLA",
        options={"maxiter": 200, "disp": True},
    )
session.close()
print(result)

Крок 4: Постобробка, повернення результату в класичному форматі

Ми бачимо, що процедура мінімізації успішно завершилась — тобто ми досягли допуску за замовчуванням класичного оптимізатора COBYLA. Якщо потрібен точніший результат, можна вказати менший допуск. Це може справді знадобитися, оскільки результат відрізнявся від отриманого симулятором вище на кілька відсотків.

Значення x, отримане в результаті, є поточним найкращим наближенням параметрів, що мінімізують функцію вартості. Якщо потрібна вища точність при подальших ітераціях, ці значення слід використовувати замість початкового x0 (вектора з одиниць).

Нарешті, зазначимо, що функція була обчислена 96 разів у процесі оптимізації. Це може відрізнятися від кількості кроків оптимізації, оскільки деякі оптимізатори вимагають кількох обчислень функції за один крок — наприклад, при оцінці градієнта.

Subspace Search VQE (SSVQE)

SSVQE — варіант VQE, що дозволяє знаходити перші $k$ власних значень спостережуваного оператора $\hat{H}$ із власними значеннями $\{\lambda_0, \lambda_1,...,\lambda_{N-1}\}$, де $N\geq k$, а не лише перше. Без обмеження загальності припускаємо, що $\lambda_0<\lambda_1<...<\lambda_{N-1}$. SSVQE вводить нову ідею: додавання ваг, щоб пріоритизувати оптимізацію за доданком із найбільшою вагою.

Для реалізації цього алгоритму потрібні $k$ взаємно ортогональних референсних станів $\{ |\rho_j\rangle \}_{j=0}^{k-1}$, тобто $\langle \rho_j | \rho_l \rangle = \delta_{jl}$ для $j,l<k$. Ці стани можна побудувати за допомогою операторів Паулі. Функція вартості цього алгоритму тоді матиме вигляд:

$$\begin{aligned} C(\vec{\theta}) & := \sum_{j=0}^{k-1} w_j \langle \rho_j | U_{V}^{\dagger}(\vec{\theta})\hat{H} U_{V}(\vec{\theta})|\rho_j \rangle \\[1mm] & := \sum_{j=0}^{k-1} w_j \langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle \\[1mm] \end{aligned}$$

де $w_j$ — довільне додатне число таке, що якщо $j<l<k$, то $w_j>w_l$, а $U_V(\vec{\theta})$ — задана користувачем варіаційна форма.

Алгоритм SSVQE спирається на той факт, що власні стани, що відповідають різним власним значенням, взаємно ортогональні. Зокрема, скалярний добуток $U_V(\vec{\theta})|\rho_j\rangle$ та $U_V(\vec{\theta})|\rho_l\rangle$ можна записати як:

$$\begin{aligned} \langle \rho_j | U_{V}^{\dagger}(\vec{\theta})U_{V}(\vec{\theta})|\rho_l \rangle & = \langle \rho_j | I |\rho_l \rangle \\[1mm] & = \langle \rho_j | \rho_l \rangle \\[1mm] & = \delta_{jl} \end{aligned}$$

Перша рівність виконується, оскільки $U_{V}(\vec{\theta})$ — квантовий оператор і тому є унітарним. Остання рівність виконується через ортогональність референсних станів $|\rho_j\rangle$. Той факт, що ортогональність зберігається при унітарних перетвореннях, глибоко пов'язаний із принципом збереження інформації, як його формулює квантова інформатика. З цієї точки зору, неунітарні перетворення відповідають процесам, де інформація втрачається або додається.

Ваги $w_j$ допомагають забезпечити, щоб усі стани ставали власними. Якщо ваги достатньо різні, доданок із найбільшою вагою (тобто $w_0$) отримуватиме пріоритет під час оптимізації над іншими. В результаті стан $U_{V}(\vec{\theta})|\rho_0 \rangle$ стає власним станом, що відповідає $\lambda_0$. Оскільки $\{ U_{V}(\vec{\theta})|\rho_j\rangle \}_{j=0}^{k-1}$ взаємно ортогональні, решта станів будуть ортогональними до нього і, відповідно, міститимуться в підпросторі, що відповідає власним значенням $\{\lambda_1,...,\lambda_{N-1}\}$.

Застосовуючи той самий аргумент до решти доданків, наступним пріоритетом буде доданок із вагою $w_1$, тому $U_{V}(\vec{\theta})|\rho_1 \rangle$ стане власним станом для $\lambda_1$, а решта доданків міститимуться у власному просторі $\{\lambda_2,...,\lambda_{N-1}\}$.

Міркуючи індуктивно, ми доходимо висновку, що $U_{V}(\vec{\theta})|\rho_j \rangle$ буде наближеним власним станом для $\lambda_j$ при $0\leq j < k.$

Теоретична схема

Схему SSVQE можна підсумувати так:

• Підготуй кілька референсних станів, застосовуючи унітарний оператор U_R до $k$ різних обчислювальних базисних станів
  • Цей алгоритм вимагає використання $k$ взаємно ортогональних референсних станів $\{ |\rho_j\rangle \}_{j=0}^{k-1}$, таких що $\langle \rho_j | \rho_l \rangle = \delta_{jl}$ для $j,l<k$.
• Застосуй варіаційну форму $U_V(\vec\theta_{i,j})$ до кожного референсного стану, отримуючи такий анзац $U_A(\vec\theta_{i,j})$.
• Виконай ініціалізацію при $i=0$, якщо доступна подібна задача (зазвичай знайдена через класичне моделювання або вибірку).
• Обчисли функцію вартості $C(\vec\theta_{i,j}) := \sum_{j=0}^{k-1} w_j \langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ для всіх підготовлених станів на квантовому комп'ютері.
  • Це можна розділити на обчислення математичного сподівання для спостережуваного оператора $\langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ та множення результату на $w_j$.
  • Після цього функція вартості повертає суму всіх зважених математичних сподівань.
• Використай класичний оптимізатор для визначення наступного набору параметрів $\Theta_{i+1}$.
• Повторюй наведені кроки до досягнення збіжності.

Ти будеш відтворювати функцію вартості SSVQE в завданні для оцінювання, але ось фрагмент для натхнення:

import numpy as np

def cost_func_ssvqe(
    params, initialized_anastz_list, weights, ansatz, hamiltonian, estimator
):
    # """Return estimate of energy from estimator

    # Parameters:
    #     params (ndarray): Array of ansatz parameters
    #     initialized_anastz_list (list QuantumCircuit): Array of initialised ansatz with reference
    #     weights (list): List of weights
    #     ansatz (QuantumCircuit): Parameterized ansatz circuit
    #     hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
    #     estimator (Estimator): Estimator primitive instance

    # Returns:
    #     float: Weighted energy estimate
    # """

    energies = []

    # Define SSVQE

    weighted_energy_sum = np.dot(energies, weights)
    return weighted_energy_sum

Variational Quantum Deflation (VQD)

VQD — ітераційний метод, що розширює VQE для отримання перших $k$ власних значень спостережуваного оператора $\hat{H}$ із власними значеннями $\{\lambda_0, \lambda_1,...,\lambda_{N-1}\}$, де $N\geq k$, а не лише першого. Далі в цьому розділі ми без обмеження загальності припускаємо, що $\lambda_0\leq\lambda_1\leq...\leq\lambda_{N-1}$. VQD вводить поняття штрафної вартості для керування процесом оптимізації.

VQD вводить штрафний доданок, позначений як $\beta$, для балансування внеску кожного доданка перекриття у функцію вартості. Цей штраф карає процес оптимізації, якщо ортогональність не досягається. Ми накладаємо це обмеження, бо власні стани спостережуваного оператора (або ермітового оператора), що відповідають різним власним значенням, завжди взаємно ортогональні або можуть бути зроблені такими у випадку виродження чи кратних власних значень. Таким чином, примусово забезпечуючи ортогональність до власного стану для $\lambda_0$, ми фактично оптимізуємо у підпросторі, що відповідає решті власних значень $\{\lambda_1, \lambda_2,..., \lambda_{N-1}\}$. При цьому $\lambda_1$ є найменшим власним значенням серед решти, і тому оптимальне розв'язання нової задачі можна отримати за допомогою варіаційної теореми.

Загальна ідея VQD полягає у використанні VQE у звичайний спосіб для знаходження найменшого власного значення $\lambda_0 := C_0(\vec\theta^0) \equiv C_\text{VQE}(\vec\theta^0)$ разом із відповідним (наближеним) власним станом $|\psi(\vec{\theta^0})\rangle$ для деякого оптимального вектора параметрів $\vec{\theta^0}$. Потім, щоб знайти наступне власне значення $\lambda_1 > \lambda_0$, замість мінімізації функції вартості $C_0(\vec{\theta}) := \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$, ми оптимізуємо:

$$C_1(\vec{\theta}) := C_0(\vec{\theta})+ \beta_0 |\langle \psi(\vec{\theta})| \psi(\vec{\theta^0})\rangle |^2$$

Додатне значення $\beta_0$ в ідеалі має бути більшим за $\lambda_1-\lambda_0$.

Це вводить нову функцію вартості, яку можна розглядати як задачу з обмеженнями: ми мінімізуємо $C_\text{VQE}(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ за умови, що стан має бути ортогональним до раніше знайденого $|\psi(\vec{\theta^0})\rangle$, а $\beta_0$ є штрафом за невиконання цього обмеження.

Альтернативно, цю нову задачу можна інтерпретувати як запуск VQE для нового спостережуваного оператора:

$$\hat{H_1} := \hat{H} + \beta_0 |\psi(\vec{\theta^0})\rangle \langle \psi(\vec{\theta^0})| \quad \Rightarrow \quad C_1(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H_1} | \psi(\vec{\theta})\rangle,$$

Припускаючи, що розв'язком нової задачі є $|\psi(\vec{\theta^1})\rangle$, математичне сподівання $\hat{H}$ (не $\hat{H_1}$) повинно дорівнювати $\langle \psi(\vec{\theta^1}) | \hat{H} | \psi(\vec{\theta^1})\rangle = \lambda_1$. Для отримання третього власного значення $\lambda_2$ функція вартості, яку треба оптимізувати, має вигляд:

$$C_2(\vec{\theta}) := C_1(\vec{\theta}) + \beta_1 |\langle \psi(\vec{\theta})| \psi(\vec{\theta^1})\rangle |^2$$

де $\beta_1$ — достатньо велика додатна константа для забезпечення ортогональності розв'язкового стану до обох $|\psi(\vec{\theta^0})\rangle$ та $|\psi(\vec{\theta^1})\rangle$. Це штрафує стани в просторі пошуку, які не відповідають цій вимозі, фактично обмежуючи простір пошуку. Таким чином, оптимальним розв'язком нової задачі повинен бути власний стан для $\lambda_2$.

Як і попередній випадок, цю нову задачу також можна інтерпретувати як VQE зі спостережуваним оператором:

$$\hat{H_2} := \hat{H_1} + \beta_1 |\psi(\vec{\theta^1})\rangle \langle \psi(\vec{\theta^1})| \quad \Rightarrow \quad C_2(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H_2} | \psi(\vec{\theta})\rangle.$$

Якщо розв'язком цієї нової задачі є $|\psi(\vec{\theta^2})\rangle$, математичне сподівання $\hat{H}$ (не $\hat{H_2}$) повинно дорівнювати $\langle \psi(\vec{\theta^2}) | \hat{H} | \psi(\vec{\theta^2})\rangle = \lambda_2$. Аналогічно, для отримання $k$-го власного значення $\lambda_{k-1}$ треба мінімізувати функцію вартості:

$$C_{k-1}(\vec{\theta}) := C_{k-2}(\vec{\theta}) + \beta_{k-2} |\langle \psi(\vec{\theta})| \psi(\vec{\theta^{k-2}})\rangle |^2,$$

Нагадаємо, що ми визначили $\vec{\theta^j}$ так, що $\langle \psi(\vec{\theta^j}) | \hat{H} | \psi(\vec{\theta^j})\rangle = \lambda_j, \forall j<k$. Ця задача еквівалентна мінімізації $C(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ за умови, що стан має бути ортогональним до $|\psi(\vec{\theta^j})\rangle ; \forall j \in {0, \cdots, k-1}$, тобто простір пошуку обмежений підпростором, що відповідає власним значенням $\{\lambda_{k-1},\cdots,\lambda_{N-1}\}$.

Ця задача еквівалентна VQE зі спостережуваним оператором:

$$\hat{H}_{k-1} := \hat{H}_{k-2} + \beta_{k-2} |\psi(\vec{\theta^{k-2}})\rangle \langle \psi(\vec{\theta^{k-2}})| \quad \Rightarrow \quad C_{k-1}(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H}_{k-1} | \psi(\vec{\theta})\rangle,$$

Як бачиш із цього процесу, для отримання $k$-го власного значення потрібні (наближені) власні стани попередніх $k-1$ власних значень, тому VQE доведеться запустити загалом $k$ разів. Отже, функція вартості VQD має вигляд:

$$C_k(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle + \sum_{j=0}^{k-1}\beta_j |\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2$$

Теоретична схема

Схему VQD можна підсумувати так:

• Підготуй референсний оператор $U_R$
• Застосуй варіаційну форму $U_V(\vec\theta_{i,j})$ до референсного стану, створюючи анзаци $U_A(\vec\theta_{i,j})$
• Виконай ініціалізацію при $i=0$, якщо є подібна задача (зазвичай знайдена через класичне моделювання або вибірку).
• Обчисли функцію вартості $C_k(\vec{\theta})$, що включає розрахунок $k$ збуджених станів та масиву значень $\beta$, що визначають штраф перекриття для кожного доданка.
  • Обчисли математичне сподівання для спостережуваного оператора $\langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ для кожного $k$
  • Обчисли штраф $\sum_{j=0}^{k-1}\beta_j |\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2$.
  • Функція вартості повертає суму цих двох доданків
• Використай класичний оптимізатор для вибору наступного набору параметрів $\Theta_{i+1}$.
• Повторюй цей процес до досягнення збіжності.

Реалізація

Для цієї реалізації ми створимо функцію для штрафу перекриття. Цей штраф буде використовуватися в функції вартості на кожній ітерації. Цей процес повторюватиметься для кожного збудженого стану.

from qiskit.circuit.library import TwoLocal

ansatz = TwoLocal(2, rotation_blocks=["ry", "rz"], entanglement_blocks="cz", reps=1)

ansatz.decompose().draw("mpl")

Спочатку налаштуємо функцію, що обчислює вірність стану — відсоток перекриття між двома станами, який використовуватиметься як штраф для VQD:

import numpy as np

def calculate_overlaps(ansatz, prev_circuits, parameters, sampler):
    def create_fidelity_circuit(circuit_1, circuit_2):
        """
        Constructs the list of fidelity circuits to be evaluated.
        These circuits represent the state overlap between pairs of input circuits,
        and their construction depends on the fidelity method implementations.
        """

        if len(circuit_1.clbits) > 0:
            circuit_1.remove_final_measurements()
        if len(circuit_2.clbits) > 0:
            circuit_2.remove_final_measurements()

        circuit = circuit_1.compose(circuit_2.inverse())
        circuit.measure_all()
        return circuit

    overlaps = []

    for prev_circuit in prev_circuits:
        fidelity_circuit = create_fidelity_circuit(ansatz, prev_circuit)
        sampler_job = sampler.run([(fidelity_circuit, parameters)])
        meas_data = sampler_job.result()[0].data.meas

        counts_0 = meas_data.get_int_counts().get(0, 0)
        shots = meas_data.num_shots
        overlap = counts_0 / shots
        overlaps.append(overlap)

    return np.array(overlaps)

Тепер час написати функцію вартості для VQD. Як і раніше при обчисленні лише основного стану, ми визначимо стан із найменшою енергією за допомогою примітиву Estimator. Однак, як описано вище, тепер ми додамо штрафний доданок для забезпечення ортогональності станів із вищою енергією. Тобто для кожного нового збудженого стану додається штраф за будь-яке перекриття між поточним варіаційним станом та вже знайденими власними станами з нижчою енергією.

def cost_func_vqd(
    parameters, ansatz, prev_states, step, betas, estimator, sampler, hamiltonian
):
    estimator_job = estimator.run([(ansatz, hamiltonian, [parameters])])

    total_cost = 0

    if step > 1:
        overlaps = calculate_overlaps(ansatz, prev_states, parameters, sampler)
        total_cost = np.sum(
            [np.real(betas[state] * overlap) for state, overlap in enumerate(overlaps)]
        )

    estimator_result = estimator_job.result()[0]

    value = estimator_result.data.evs[0] + total_cost

    return value

Зверни увагу, що наведена вище функція вартості звертається до функції calculate_overlaps, яка фактично створює новий квантовий Circuit. Якщо ми хочемо запустити це на реальному залізі, цей новий Circuit також має бути транспільований — бажано оптимальним чином — для Backend, який ми обираємо. Зауваж, що транспіляцію не вбудовано у функції calculate_overlaps або cost_func_vqd. Спробуй самостійно модифікувати код, щоб додати цю додаткову (і умовну) транспіляцію — але в наступному уроці це буде зроблено за тебе.

У цьому уроці ми запустимо алгоритм VQD за допомогою Statevector Sampler та Statevector Estimator:

from qiskit.primitives import StatevectorEstimator as Estimator

sampler = Sampler()
estimator = Estimator()

Введемо спостережуваний оператор для оцінювання. У наступному уроці ми додамо до цього фізичний контекст — наприклад, збуджений стан молекули. Корисно уявляти цей спостережуваний оператор як гамільтоніан системи, що може мати збуджені стани, хоча його не підібрано для відповідності якійсь конкретній молекулі чи атому.

from qiskit.quantum_info import SparsePauliOp

observable = SparsePauliOp.from_list([("II", 2), ("XX", -2), ("YY", 3), ("ZZ", -3)])

Тут ми задаємо загальну кількість станів, які хочемо обчислити (основний стан та збуджені стани, k), і штрафи (betas) за перекриття між векторами стану, які мають бути ортогональними. Наслідки вибору надто великих або надто малих betas будуть досліджені в наступному уроці. Наразі ми просто використаємо наведені нижче значення. Починатимемо з усіх нулів як параметрів. У власних обчисленнях ти можеш захотіти використовувати більш розумні початкові параметри, спираючись на знання системи або попередні розрахунки.

k = 3
betas = [33, 33, 33]
x0 = np.zeros(8)

Тепер можна запустити обчислення:

from scipy.optimize import minimize

prev_states = []
prev_opt_parameters = []
eigenvalues = []

for step in range(1, k + 1):
    if step > 1:
        prev_states.append(ansatz.assign_parameters(prev_opt_parameters))

    result = minimize(
        cost_func_vqd,
        x0,
        args=(ansatz, prev_states, step, betas, estimator, sampler, observable),
        method="COBYLA",
        options={
            "maxiter": 200,
        },
    )
    print(result)

    prev_opt_parameters = result.x
    eigenvalues.append(result.fun)

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -5.999999979545955
       x: [-5.150e-01 -5.452e-02 -1.571e+00 -2.853e-05  2.671e-01
           -2.672e-01 -8.509e-01 -8.510e-01]
    nfev: 131
   maxcv: 0.0
 message: Optimization terminated successfully.
 success: True
  status: 1
     fun: 4.024550284767612
       x: [-3.745e-01  1.041e+00  8.637e-01  1.202e+00 -8.847e-02
            1.181e-02  7.611e-01 -3.006e-01]
    nfev: 110
   maxcv: 0.0
 message: Optimization terminated successfully.
 success: True
  status: 1
     fun: 5.608925562838559
       x: [-2.670e-01  1.280e+00  1.070e+00 -8.031e-01 -1.524e-01
           -6.956e-02  7.018e-01  1.514e+00]
    nfev: 90
   maxcv: 0.0

Отримані значення функції вартості становлять приблизно -6.00, 4.02 і 5.61. Найважливіше в цих результатах — те, що значення функції зростають. Якби ми отримали перший збуджений стан із нижчою енергією, ніж у початковому необмеженому розрахунку основного стану, це вказувало б на помилку в коді.

Значення x — це параметри, що дали вектор стану, який відповідає кожній із цих вартостей (енергій).

Нарешті, зазначимо, що всі три мінімізації збіглися в межах допуску за замовчуванням класичного оптимізатора (у даному випадку COBYLA). Для них знадобилося 131, 110 і 90 обчислень функції відповідно.

Квантова регресія з вибіркою (QSR)

Одна з головних проблем VQE — це численні звернення до квантового комп'ютера, необхідні для отримання параметрів на кожному кроці: наприклад, $k$, $k-1$ і так далі. Це особливо незручно, коли доступ до квантових пристроїв ставиться в чергу. Хоча Session дозволяє об'єднати кілька ітеративних викликів, альтернативним підходом є використання вибірки. Задіявши більше класичних ресурсів, можна виконати весь процес оптимізації в одному зверненні. Саме тут у гру вступає Квантова регресія з вибіркою (Quantum Sampling Regression). Оскільки доступ до квантових комп'ютерів досі залишається ресурсом із низькою пропозицією та високим попитом, такий компроміс видається водночас можливим і зручним для багатьох сучасних досліджень. Цей підхід максимально використовує всі наявні класичні обчислювальні можливості й водночас зберігає чимало внутрішніх механізмів та власних властивостей квантових обчислень, які не відтворюються в симуляції.

Ідея QSR полягає в тому, що функція витрат $C(\theta) := \langle \psi(\theta) | \hat{H} | \psi(\theta)\rangle$ може бути представлена у вигляді ряду Фур'є таким чином:

$$\begin{aligned} C(\vec{\theta}) & := \langle \psi(\theta) | \hat{H} | \psi(\theta)\rangle \\[1mm] & := a_0 + \sum_{k=1}^S[a_k\cos(k\theta)+ b_k\sin(k\theta)] \\[1mm] \end{aligned}$$

Залежно від periodicity та смуги пропускання вихідної функції, множина $S$ може бути скінченною або нескінченною. Для цілей нашого обговорення припускатимемо, що вона нескінченна. Наступний крок — кількаразово відібрати вибірки функції витрат $C(\theta)$, щоб отримати коефіцієнти Фур'є $\{a_0, a_k, b_k\}_{k=1}^S$. Зокрема, оскільки ми маємо $2S+1$ невідомих, нам потрібно відібрати вибірки функції витрат рівно $2S+1$ разів.

Якщо ми відберемо вибірки функції витрат для $2S+1$ значень параметра $\{\theta_1,...,\theta_{2S+1}\}$, отримаємо таку систему:

$$\begin{pmatrix} 1 & \cos(\theta_1) & \sin(\theta_1) & \cos(2\theta_1) & ... & \sin(S\theta_1) \\ 1 & \cos(\theta_2) & \sin(\theta_2) & \cos(2\theta_2) & \cdots & \sin(S\theta_2)\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & \cos(\theta_{2S+1}) & \sin(\theta_{2S+1}) & \cos(2\theta_{2S+1}) & \cdots & \sin(S\theta_{2S+1}) \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ b_1 \\ a_2 \\ \vdots \\ b_S \end{pmatrix} = \begin{pmatrix} C(\theta_1) \\ C(\theta_2) \\ \vdots \\ C(\theta_{2S+1}) \end{pmatrix},$$

яку перепишемо як

$$Fa=c.$$

На практиці ця система, як правило, не є сумісною, оскільки значення функції витрат $c$ не є точними. Тому зазвичай варто нормалізувати їх, помноживши зліва на $F^\dagger$, що дає:

$$F^\dagger Fa = F^\dagger c.$$

Ця нова система завжди є сумісною, а її розв'язок є розв'язком задачі найменших квадратів для вихідної проблеми. Якщо ми маємо $k$ параметрів замість одного, і кожен параметр $\theta^i$ має власне $S_i$ для $i \in {1,...,k}$, то загальна кількість необхідних вибірок становить:

$$T=\prod_{i=1}^k(2S_i+1)\leq \prod_{i=1}^k(2S_{max}+1) = (2S_{max}+1)^n,$$

де $S_{\max} = \max_i(S_i)$. Крім того, якщо задавати $S_{\max}$ як настроюваний параметр (а не виводити його), відкриваються нові можливості, зокрема:

• Надмірна вибірка (over-sampling) для підвищення точності.
• Недостатня вибірка (under-sampling) для прискорення роботи шляхом зменшення накладних витрат або усунення локальних мінімумів.

Теоретична схема

Схему QSR можна описати так:

• Підготувати референсні оператори $U_R$.
  • Ми переходимо зі стану $|0\rangle$ до референсного стану $|\rho\rangle$
• Застосувати варіаційну форму $U_V(\vec\theta_{i,j})$, щоб створити анзац $U_A(\vec\theta_{i,j})$.
  • Визначити смугу пропускання для кожного параметра в анзаці. Достатньо верхньої межі.
• Виконати bootstrap при $i=0$, якщо є подібна задача (як правило, знайдена класичною симуляцією або вибіркою).
• Відібрати вибірки функції витрат $C(\vec\theta) := a_0 + \sum_{k=1}^S[a_k\cos(k\theta)+ b_k\sin(k\theta)]$ щонайменше $T$ разів.
  • $T=\prod_{i=1}^k(2S_i+1)\leq \prod_{i=1}^k(2S_{max}+1) = (2S_{max}+1)^n$
  • Вирішити, чи використовувати надмірну/недостатню вибірку, щоб збалансувати швидкість і точність, регулюючи $T$.
• Обчислити коефіцієнти Фур'є з вибірок (тобто розв'язати нормалізовану лінійну систему рівнянь).
• Знайти глобальний мінімум отриманої функції регресії на класичній машині.

Підсумок

З цього уроку ти дізнався про кілька варіаційних підходів:

• Загальна схема
• Введення ваг і штрафів для коригування функції витрат
• Дослідження недостатньої та надмірної вибірки для компромісу між швидкістю й точністю

Ці ідеї можна адаптувати для створення власного варіаційного алгоритму, що відповідає твоїй задачі. Заохочуємо тебе ділитися своїми результатами зі спільнотою. У наступному уроці ми розглянемо, як використовувати варіаційний алгоритм для розв'язання прикладної задачі.