Практика VQE
Contents
Практика VQE#
Автор(ы):
Эта лекция будет посвящена практике с VQE, а именно решению задачи коммивояжера. Будет много кода и мало теории. Рекомендуется заранее повторить следующие лекции:
Описание задачи коммивояжера в лекции про задачи комбинаторной оптимизации.
Сведение этой задачи к модели Изинга в лекции о переходе от комбинаторных задач к квантовым алгоритмам.
Нам потребуется определить бинарный оператор
Импорт библиотек#
Для графов будем использовать структуры данных из уже знакомой нам библиотеки NetwrokX
, а для работы с разреженными операторами SciPy
.
import matplotlib.pyplot as plt
%config InlineBackend.figure_format = "retina"
import numpy as np
import networkx as nx
from scipy import sparse
from scipy.sparse import linalg as sl
Случайный взвешенный граф#
Случайный взвешенный граф проще всего сгенерировать из случайной матрицы весов. Давайте сделаем такую:
Note
Тут будем работать с матрицей полного графа, то есть такого, в котором каждая вершина связана с каждой. Это вполне корректное допущение для задачи коммивояжера, так как, например, из каждого города РФ обычно можно попасть во все другие, пусть и за разное время в пути.
np.random.seed(42)
rand_mat = np.random.randint(0, 100, (4, 4))
Только эта матрица не является симметричной и имеет ненулевые элементы на диагонали. Элементы на диагонали означают петли (loop или self-loop), что не очень разумно в терминах задачи коммивояжера: не ясно, что значит, например, путь из Москвы в Москву. Аналогично и про то, что матрица не симметричная: если путь из Москвы в Санкт-Петербург займет 6 часов, то и путь обратно тоже займет столько же. Так что хотим симметричную матрицу. Давайте это исправим:
rand_adj = (rand_mat + rand_mat.T) / 2
np.fill_diagonal(rand_adj, 0)
Наконец, сделаем из этого граф и визуализируем его:
g = nx.Graph(rand_adj)
pos = nx.drawing.layout.random_layout(g, seed=42)
nx.draw(g, pos=pos)
edge_weights = nx.get_edge_attributes(g, "weight")
nx.draw_networkx_edge_labels(g, pos=pos, edge_labels=edge_weights)
nx.draw_networkx_labels(g, pos=pos)
plt.show()

Решение задачи TSP в NetworkX#
Для того, чтобы потом понять, что все правильно сделали, давайте найдем приближенное решение задачи средствами NetworkX
.
nx_solution = nx.algorithms.approximation.traveling_salesman_problem(g)
print(f"NetworkX solution (nodes order): {nx_solution}")
nx_solution_cost = 0
for i in range(len(nx_solution) - 1):
nx_solution_cost += g.get_edge_data(nx_solution[i], nx_solution[i + 1])["weight"]
print(f"NetworkX solution (total cost): {nx_solution_cost}")
NetworkX solution (nodes order): [0, 3, 1, 2, 0]
NetworkX solution (total cost): 213.0
Матрица #
Это бинарный оператор, который ранее вводили в лекции о представлении комбинаторных задач в виде моделей Изинга.
Соотношение между
def sigma_z_mod(k: int, size: int) -> sparse.csr_matrix:
sz_mod = sparse.coo_matrix(
[
[1, 0,],
[0, 0,],
]
)
left_ith = sparse.eye(2 ** k, format="coo")
right_ith = sparse.eye(2 ** (size - k - 1), format="coo")
return sparse.kron(
sparse.kron(left_ith, sz_mod),
right_ith,
)
Для произвольного состояния
Пусть
Отсюда следует, что
Значит,
собственные векторы
и совпадают;соответствующие собственные значения связаны соотношением
.
Поскольку
eigs = np.linalg.eig(sigma_z_mod(1, 3).toarray())
print(np.unique(eigs[0]))
[0. 1.]
Напомню, что у нас zero-based индексация, поэтому sigma_z_mod(1, 3)
значит что-то типа “построй мне матрицу
Легко убедиться, что
для
собственный вектор представляется битовой строкой, в которой спин с индексом 1 имеет значение 0;для
собственный вектор представляется битовой строкой, в которой спин с индексом 1 имеет значение 1.
def sigmaz_k(k: int, n: int) -> sparse.csr_matrix:
left_part = sparse.eye(2 ** k)
right_part = sparse.eye(2 ** (n - 1 - k))
return sparse.kron(
sparse.kron(
left_part,
sparse.csr_matrix(np.array([[1, 0,], [0, -1,],]))
),
right_part
)
def probs2bit_str(probs: np.array) -> str:
size = int(np.log2(probs.shape[0]))
bit_s_num = np.where(probs == probs.max())[0][0]
s = f"{bit_s_num:b}"
s = "0" * (size - len(s)) + s
return s
for i, psi in enumerate(eigs[1]):
probs = psi * psi.conj()
bit_s = probs2bit_str(probs)
print(f"Eigenvalue: {eigs[0][i]}\tBit string: {bit_s}")
Eigenvalue: 1.0 Bit string: 000
Eigenvalue: 1.0 Bit string: 001
Eigenvalue: 0.0 Bit string: 010
Eigenvalue: 0.0 Bit string: 011
Eigenvalue: 1.0 Bit string: 100
Eigenvalue: 1.0 Bit string: 101
Eigenvalue: 0.0 Bit string: 110
Eigenvalue: 0.0 Bit string: 111
Убедимся также в совпадении собственных векторов с матрицей
Note
Для чистых состояний без запутывания матрицы собственных векторов могут быть приведены к виду единичных матриц. И действительно, каждый собственный вектор отвечает строго одной конфигурации спинов, а значит его волновая функция имеет вид
print(f"Sigma z_mod:\n{eigs[1]}")
print(f"\nSigma z:\n{np.linalg.eig(sigmaz_k(1, 3).toarray())[1]}")
Sigma z_mod:
[[1. 0. 0. 0. 0. 0. 0. 0.]
[0. 1. 0. 0. 0. 0. 0. 0.]
[0. 0. 1. 0. 0. 0. 0. 0.]
[0. 0. 0. 1. 0. 0. 0. 0.]
[0. 0. 0. 0. 1. 0. 0. 0.]
[0. 0. 0. 0. 0. 1. 0. 0.]
[0. 0. 0. 0. 0. 0. 1. 0.]
[0. 0. 0. 0. 0. 0. 0. 1.]]
Sigma z:
[[1. 0. 0. 0. 0. 0. 0. 0.]
[0. 1. 0. 0. 0. 0. 0. 0.]
[0. 0. 1. 0. 0. 0. 0. 0.]
[0. 0. 0. 1. 0. 0. 0. 0.]
[0. 0. 0. 0. 1. 0. 0. 0.]
[0. 0. 0. 0. 0. 1. 0. 0.]
[0. 0. 0. 0. 0. 0. 1. 0.]
[0. 0. 0. 0. 0. 0. 0. 1.]]
Note
В общем случае, даже если у двух операторов совпадают собственные вектора, то их порядок (по возрастанию собственных значений) совпадать совсем не обязан. В этом легко убедиться, например, если поменять местами np.allclose
, а сравнивать их глазами построчно.
Гамильтониан задачи TSP#
Давайте еще раз выпишем полный гамильтониан системы. Напомню, что задача будет параметризована при помощи
Note
Еще раз вспомним, что первые два слагаемых “говорят”: каждая вершина должна входить только в один столбец, а каждая строка должна содержать только одну вершину. Третье слагаемое представляет собой штраф за соединение вершин, между которыми нет ребра, а четвертое дает стоимость всего пути. Первые три слагаемых идут с бОльшим коэффициентом, чем последнее, так как нам важно в первую очередь получить реальное решение, а не получить самое дешевое.
Сделаем финт ушами: будем всегда считать, что путь коммивояжера начинается в вершине с индексом 0. Так как путь коммивояжера это гамильтонов цикл, то делаем это упрощение без какой-либо потери общности, так как вершина с индексом 0 все равно должна будет войти в путь. При этом такой хитрый трюк позволит снизить число необходимых кубит с NetworkX
:
def convert_tsp_to_ising(g: nx.Graph) -> sparse.coo_matrix:
# Имеется два элемента гамильтониана:
# - часть, дающая нам Гамилтонов цикл
# - часть, считающая стоимость пути
# Для нас важно, чтобы вес первой части превалирвал.
# Поэтому, в качестве коэффициента берем
# пятикратное значение максимального веса ребер в графе.
a = 5 * max([e[2]["weight"] for e in g.edges(data=True)])
size = g.number_of_nodes() - 1
H = sparse.coo_matrix((2 ** (size ** 2), 2 ** (size ** 2)), dtype=np.complex128)
# Вычисляем сумму по колонкам:
# - каждая вершина должна быть в каждой колонке строго один раз
for i in range(size):
# Считаем выражение 1 - sum
col_sum = sparse.eye(2 ** (size ** 2), format="coo")
for j in range(size):
col_sum -= sigma_z_mod(i + j * size, size ** 2)
# Возводим в квадрат и добавляем к гамилтониану
H += a * col_sum * col_sum
# Вычисляем сумму по строкам:
# - каждая вершина должна быть в каждой строке строго один раз
for j in range(size):
row_sum = sparse.eye(2 ** (size ** 2), format="coo")
for i in range(size):
row_sum -= sigma_z_mod(i + j * size, size ** 2)
# Аналогично, в квадрат и добавляем
H += a * row_sum * row_sum
# Считаем сумму по тем ребрам, которых в графе нет и тем, которые есть.
for i in range(size):
for j in range(size):
if g.has_edge(i + 1, j + 1):
# Сумма по всем ребрам графа
coef = g.get_edge_data(i + 1, j + 1)["weight"]
else:
# Сумма по всем отсутствующим ребрам графа
coef = a
for k in range(size - 1):
# Дальше просто используем полученный коэффициент и считаем сумму
H += (
coef
* sigma_z_mod(i + k * size, size ** 2)
* sigma_z_mod(j + (k + 1) * size, size ** 2)
)
H += (
coef
* sigma_z_mod(i + (size - 1) * size, size ** 2)
* sigma_z_mod(j, size ** 2)
)
# Начинаем всегда из нулевой вершины, это позволяет снизить число необходимых кубитов.
# Теперь должны добавить в гамильтониан слагаемое, отвечающее за нулевую вершину в пути.
for j in range(size):
if g.has_edge(0, j + 1):
coef = g.get_edge_data(0, j + 1)["weight"]
else:
coef = a
H += coef * sigma_z_mod(j, size ** 2)
H += coef * sigma_z_mod(j + (size - 1) * size, size ** 2)
# PennyLane работает с матрицами только в координатном формате, поэтому сразу делаем приведение.
return sparse.coo_matrix(H)
Убедимся, что точное решение даст тот же путь (с той же стоимостью), что и решение из NetworkX
:
op = convert_tsp_to_ising(g)
scipy_solution = sl.eigs(op, k=1, which="SR", return_eigenvectors=True)
print(f"Ground state energy: {scipy_solution[0][0]:.4f}")
probs = scipy_solution[1] * scipy_solution[1].conj()
bit_s = probs2bit_str(probs)
# Нулевая вершина всегда первая
scipy_path = [0,]
for i in range(3):
k = 0
for j in range(3):
if bit_s[j + i * 3] == "0":
k = j + 1
scipy_path.append(k)
scipy_path.append(0)
scipy_solution_cost = 0
for i in range(len(scipy_path) - 1):
scipy_solution_cost += g.get_edge_data(scipy_path[i], scipy_path[i + 1])["weight"]
print(f"SciPy solution (nodes order): {scipy_path}")
print(f"SciPy solution (total cost): {scipy_solution_cost}")
Ground state energy: 273.0000-0.0000j
SciPy solution (nodes order): [0, 2, 1, 3, 0]
SciPy solution (total cost): 213.0
Видим, что решение через матрицу совпадает с тем, что выдал NetworkX
(с точностью до инверсии: [0, 2, 1, 3, 0]
и [0, 3, 1, 2, 0]
), как и общая стоимость пути. Энергия основного состояния отличается по причине того, что использовали оптимизацию и магию с коэффициентами. Но это не так важно: по энергиям важно лишь отношение порядка, потому что зная битовую строку всегда можем проверить стоимость пути за линейное время!
Теперь можно наконец-то, приступать к симуляции квантовых вычислений!
VQE на Pennylane#
Для начала давайте преобразуем матрицу к тому объекту, с которым умеет работать квантовый фреймворк.
Note
Pennylane
поддерживает также работу с представлением гамильтонианов в виде так называемых Паули-строк, но тут не будем еще раз определять весь гамильтониан, а просто перериспользуем разреженную матрицу.
Веса графа имеют разброс значений от 0 до 100. Из-за этого, при некоторых конфигурациях можем получать очень большое собственное значение гамильтониана, вплоть до нескольких тысяч. Это не очень хорошо для градиентного спуска методом paramter-shift
, поэтому для численной устойчивости давайте умножим все значения весов на
import pennylane as qml
numer_stability_const = 0.1
ham = qml.SparseHamiltonian(op * numer_stability_const, wires=list(range(9)))
/home/runner/work/qmlcourse/qmlcourse/.venv/lib/python3.8/site-packages/_distutils_hack/__init__.py:33: UserWarning: Setuptools is replacing distutils.
warnings.warn("Setuptools is replacing distutils.")
Далее необходимо выбрать ansatz. Современные квантовые процессоры обычно имеют вид решетки, типа такого:

Fig. 100 Схематичное изображение квантового компьютера Sycamore от компании Google [AAB+19]#
Поэтому, в идеале, хотелось бы делать либо однокубитные операции, либо двухкубитные операции, но над “соседними” кубитами. Так появилась идея hardware efficient ansatz [KMT+17]. Он заключается в том, что формируется несколько “слоев” VQC, где каждый такой слой содержит несколько операций вращений, а также попарные “запутывающие” операции (например, CZ, или гейт CNOT). Выглядит это примерно так:

Fig. 101 Иллюстрация hardware efficient ansatz (HWEA) из [KMT+17]#
Давайте реализуем такой слой:
size = len(ham.wires)
def layer(params):
for i in range(size):
qml.RY(params[i][0], wires=i)
qml.RZ(params[i][1], wires=i)
for i in range(size - 1):
qml.CZ(wires=[i, i + 1])
Здесь специально не повторяем основ параметризированных квантовых схем. Для тех, кто забыл, что это такое, рекомендуем вернуться к отдельной лекции по этой теме.
Теперь создадим устройство, которое будет квантовым симулятором, а также напишем код, который будет на этом устройстве запускаться. Такой код должен последовательно применить несколько HWEA слоев, после чего измерить ожидаемое значение гамильтониана в полученном состоянии.
device = qml.device("default.qubit", wires=size)
num_layers = 3
@qml.qnode(device, diff_method="parameter-shift")
def feed_forward(params, wires=size):
for k in range(3):
layer(params[k])
for i in range(size):
qml.RY(params[-1][i][0], wires=i)
qml.RZ(params[-1][i][1], wires=i)
return qml.expval(ham)
Теперь надо сгенерировать параметры и выбрать метод их оптимизации. В этот раз давайте воспользуемся алгоритмом Momentum
, так называемым “ускоренным” градиентным спуском:
opt = qml.MomentumOptimizer()
params = tuple(
tuple(
(
np.random.rand() * 2 * np.pi,
np.random.rand() * 2 * np.pi,
)
for _ in range(size)
)
for _ in range(num_layers + 1)
)
Ну и наконец можно приступать к оптимизации параметров. Все примерно также, как в базовой лекции про градиенты. То есть ничего нового тут по сути не делается.
Симуляция#
Warning
Осторожно, следующий код работает несколько минут, так как обучаем довольно тяжелую схему на 9-и кубитах! Кстати, при желании можно добиться гораздо больше точности по энергии относительно точного решения. Для этого достаточно установить число эпох
energy = [
feed_forward(params),
]
for epoch in range(85):
params, e = opt.step_and_cost(feed_forward, params)
energy.append(e)
if epoch % 5 == 0:
print(f"epoch: {epoch}\tenergy: {e:.6f}")
/home/runner/work/qmlcourse/qmlcourse/.venv/lib/python3.8/site-packages/pennylane/_grad.py:95: UserWarning: Starting with PennyLane v0.21.0, when using Autograd, inputs have to explicitly specify requires_grad=True (or the argnum argument must be passed) in order for trainable parameters to be identified.
warnings.warn(
epoch: 0 energy: 360.638506
epoch: 5 energy: 163.172963
epoch: 10 energy: 96.375950
epoch: 15 energy: 56.524352
epoch: 20 energy: 43.497589
epoch: 25 energy: 35.938162
epoch: 30 energy: 31.882151
epoch: 35 energy: 29.457526
epoch: 40 energy: 29.232270
epoch: 45 energy: 28.058306
epoch: 50 energy: 27.814023
epoch: 55 energy: 27.495886
epoch: 60 energy: 27.479515
epoch: 65 energy: 27.411348
epoch: 70 energy: 27.348231
epoch: 75 energy: 27.329692
epoch: 80 energy: 27.324815
plt.figure(figsize=(10, 8))
plt.plot(np.arange(len(energy)), energy, ".-", label="Energy by epochs")
plt.hlines(
np.real(scipy_solution[0][0]) * numer_stability_const,
0,
len(energy),
color="red",
label="Exact solution",
)
plt.xlabel("Epoch")
plt.ylabel("Energy")
plt.legend()
plt.show()

Видим, что VQE отлично сходится к основному состоянию, которое знаем из SciPy
и ARPACK
. Теперь давайте убедимся, что квантовый путь совпадет с тем, что получили из NetworkX
. Для этого нужно получить распределение вероятностей битовых строк. На реальном компьютере это можно сделать просто измеряя много раз, но на симуляторе есть возможность получить его сразу одной командой.
@qml.qnode(device, diff_method="parameter-shift")
def answer(params, wires=size):
for k in range(num_layers):
layer(params[k])
for i in range(size):
qml.RY(params[-1][i][0], wires=i)
qml.RZ(params[-1][i][1], wires=i)
return qml.probs(wires=list(range(size)))
probs = answer(params)
bit_s = probs2bit_str(probs)
# Нулевая вершина всегда первая
q_path = [0,]
for i in range(3):
k = 0
for j in range(3):
if bit_s[j + i * 3] == "0":
k = j + 1
q_path.append(k)
q_path.append(0)
q_solution_cost = 0
for i in range(len(q_path) - 1):
q_solution_cost += g.get_edge_data(q_path[i], q_path[i + 1])["weight"]
print(f"Quantum solution (nodes order): {q_path}")
print(f"Quantum solution (total cost): {q_solution_cost}")
Quantum solution (nodes order): [0, 2, 1, 3, 0]
Quantum solution (total cost): 213.0
Ура, получился тот же ответ, что и в SciPy
/NetworkX
. Значит самописный VQE работает правильно.
Заключение#
В этой лекции попрактиковались в решение реальных задач при помощи алгоритма VQE.