Проблема собственных значений#
Мы с вами узнали, что задачи комбинаторной оптимизации и квантовой химии могут быть сведены к решению проблемы поиска минимального собственного значения большого эрмитова оператора – гамильтониана. Для оптимизационных задач это осуществляется при помощи сведения к QUBO-матрице и гамильтониану типа Изинга. А для электронных орбиталей из квантовой химии можно применить преобразование Жордана-Вигнера и также перейти к спиновому гамильтониану.
Теперь перед нами встает вопрос, а как же искать основное состояние этого гамильтониана? В этой лекции рассмотрим классические методы решения этой проблемы, то есть без квантовых компьютеров. Рассмотрение этих методов и их недостатков покажет то, зачем тут так нужен будет квантовый компьютер.
О проблеме (повторение)#
Эта тема обсуждалась во вводных лекциях по линейной алгебре, в части про собственные вектора и собственные значения.
или в нотации Дирака, которая используется в области квантовых вычислений:
Таким образом, собственные вектора – это такие вектора, которые при применении оператора не меняют свое направление. Например, в примере ниже собственный вектор – это ось симметрии оператора:
Итеративные алгоритмы#
В целом, задача нахождения собственных значений является очень трудной с вычислительной точки зрения, особенно для больших матриц. Для матриц размера более, чем \(3 \times 3\) в общем случае не существует алгоритма нахождения собственных значений и собственных векторов. Однако существует несколько итеративных алгоритмов. Рассмотрим лишь два из них, причем без особых деталей, так как эти алгоритмы, а также доказательство их сходимости являются достаточно сложными.
Степенной метод#
Один из самых простых для понимания алгоритмов, который, тем не менее находит интересные применения. Суть его в том, что берем некоторый случайный вектор \(\ket<\Psi>\) и начинаем последовательно действовать на него оператором \(\hat\) (другими словами умножать, на нашу матрицу), при этом нормируя:
И так повторяем до тех пор, пока изменение вектора не будет меньше, чем некоторое заданное маленькое значение \(\epsilon\) . Когда достигли этого условия, это значит что нашли первый собственный вектор, который соответствует наибольшему собственному значению. В частном случае интересных нам эрмитовых операторов, можно так же последовательно находить все собственные вектора и собственные значения.
На самом деле, сеть интернета является графом – множеством связанных между собой вершин. А любой граф можно представить в виде большой-большой, но очень разреженной матрицы, каждый элемент которой это 1 если между соответствующими вершинами есть ребро и 0, если нет. Например, элемент \(L_\) будет 1, если между вершинами \(i\) и \(j\) есть ребро.иВ 1998-м году, Ларри Пейдж и Сергей Брин нашли очень эффективный способ подсчета первого собственного вектора этой матрицы, используя именно модификацию степенного метода. Этот алгоритм получил название PageRank , причем Page это фамилия автора, а не отсылка к веб-страницам, как можно было бы подумать. Этот алгоритм лег в основу поисковика Google, который в дальнейшем вырос в транснациональную корпорацию!
Итерация Арнольди#
Это гораздо более сложный метод, который, однако, является одним из самых эффективных применительно к разреженным матрицам [Arn51] . Объяснить его легко, к сожалению, не получится, так как алгоритм требует понимания Крыловских подпространств и других концептов из области линейной алгебры разреженных систем. Но пока достаточно лишь того, что этот алгоритм имеет очень эффективную реализацию – ARPACK, написанную в середине 90-х годов на языке FORTRAN77 . Именно эта библиотека используется “под капотом” у SciPy , а также во многих других научных пакетах. Давайте посмотрим, как она работает.
Сгенерируем большую разреженную матрицу.
import numpy as np from scipy import sparse np.random.seed(42) x = np.random.random(10000) np.random.seed(42) y = np.random.random(10000) px = np.where(x > 0.2) py = np.where(y > 0.2) num_elements = max([px[0].shape[0], py[0].shape[0]]) spmat = sparse.coo_matrix( ( (np.ones(num_elements), (px[0][:num_elements], py[0][:num_elements])) ) ) print(spmat.__repr__())
' with 7957 stored elements in COOrdinate format>
Матрица размера \(10000 \times 10000\) это большая матрица и работать с ней в “плотном” (dense) представлении было бы очень трудно. Но ARPACK позволяет найти минимальное собственное значение за доли секунд, используя разреженность матрицы:
from scipy.sparse import linalg as sl max_eigval = sl.eigs(spmat, k=1, which="LR", return_eigenvectors=False)[0] min_eigval = sl.eigs(spmat, k=1, which="SR", return_eigenvectors=False)[0] print(f"Min E: min_eigval>\nMax E: max_eigval>")
Min E: (-1.1102230246251565e-16+0j) Max E: (1.0000000000000007+0j)
Для тех кто забыл, какие параметры принимает функция eigs из scipy.linalg.spare напомним, что первый параметр это разреженная матрица, k – сколько именно собственных значений хотим получить, which указывает на собственные значения:
- SM – smallest magnitude – наименьшие по модулю числа
- LM – largest magnitude – наибольшие по модулю числа
- SR – smallers real – числа с наименьшей действительной частью
- LR – largest real – числа с наибольшей действительной частью
- SI – smallest image – числа с наименьшей мнимой частью
- LI – largest image – числа с наибольшей мнимой частью
Наконец, параметр return_eigenvectors – хотим ли получить только собственные значения, или еще и собственные вектора.
Более подробна работа с scipy.sparse , а также с scipy.sparse.linalg разбирается в [вводном блоке по линейной алгебре](пока пусто).
Не у всех матриц все собственные значения являются действительными, поэтому ARPACK по умолчанию считает комплексные значения, хотя в этом конкретном случае видим, что мнимая часть равна нулю.
Алгоритм Ланкзоша#
Итерация Ланкзоша (англ. Lanzos) [Lan50] – это модификация итерации Арнольди, которая работает с эрмитовыми матрицами и находит максимально широкое применение в том числе для квантовых гамильтонианов. Этот алгоритм по умолчанию включен в большинство математических пакетов, включая ARPACK и, соответственно, SciPy :
max_eigval = sl.eigsh(spmat, k=1, which="LM", return_eigenvectors=False)[0] min_eigval = sl.eigsh(spmat, k=1, which="SM", return_eigenvectors=False)[0] print(f"Min E: min_eigval>\nMax E: max_eigval>")
Min E: -8.323011768995762e-25 Max E: 1.0000000000000002
У этой процедуры из ARPACK немного другие варианты параметра which , так как мы помним, что у эрмитовых матриц собственные значения вещественны:
- LM – largest magnitude – наибольшие по модулю
- SM – smallest magnitude – наименьшие по модулю
- LA – largest algebraic – алгебраически наибольшие, т.е. с учетом знака
- SA – smallest algebraic – алгебраически наименьшие, т.е. с учетом знака
Вариационные алгоритмы#
В этом разделе поговорим о существующих алгоритмах решения задачи об основном состоянии уже в контексте квантовой механики. Хотя, как помним, задачи оптимизации и квантовой физики тесно связаны. В каком-то смысле, вариационные алгоритмы, а в особенности, квантовый Монте-Карло и различные его модификации в чем-то сильно похожи на классический алгоритм имитации отжига.
Вариационный Монте-Карло#
Variational Monte-Carlo, или просто VMC это очень простой и в тоже время эффективный алгоритм нахождения основного состояния квантомеханической системы.
Замечание – в классическом VMC обычно работают при нулевой температуре. Хотя в общем случае, температура оказывает значительное влияние на то, в каком состоянии находится физическая система.
Давайте еще раз запишем ожидаемое значение энергии гамильтониана в состоянии \(\ket<\Psi>\) :
Если ввести вектор \(X\) , который описывает конфигурацию системы (например, ориентации спинов), то выражение для энергии можно переписать в интегральной форме:
В данном случае, выражение
дает распределение вероятностей, а значит можно из него семплировать, используя методы Монте-Карло. Это очень похоже на то, как ранее семплировали из распределения Больцмана в классическом методе Монте-Карло. Вопрос лишь в том, как представить волновую функцию \(\ket<\Psi>\) ? В этом помогут так называемые trial wave functions – параметризированные функции от \(X\) . В этом случае меняем или варьируем параметры trial wave function в процессе:
- семплируем из \(\frac<|\Psi(X)|^2><\int |\Psi(X)|^2 dX>\) конфигурации;
- обновляем параметризацию trial function так, чтобы минимизировать энергию.
Повторяем до сходимости. Ну а дальше посмотрим на некоторые примеры trial wave functions.
Jastrow Function#
Когда есть задача из \(N\) квантовых частиц, каждая из которых описывается координатой или радиус вектором, то можно построить trial wave function в виде суммы попарных функций двухчастичных взаимодействий:
где \(r_i, r_j\) – это радиус-векторы частиц, а \(u(r_i, r_j)\) – симметричная функция, описывающая двухчастичное взаимодействия. Такая функция называется Jastrow function [Jas55] . В этом случае, в процессе работы VMC будем просто варьировать радиус-векторы частиц также, как варьировали вершины графа в обычном отжиге, когда решали задачу комбинаторной оптимизации. Только теперь есть еще и параметризация обменных взаимодействий, которую “варьируем”.
Hartree-Fock (SCF)#
Для задач квантовой химии, когда работаем с фермионами, существует вид trial wave function на основе Слэтеровского детерминанта, о котором писали в продвинутой лекции по квантовой химии:
где \(D\) это матрица из одноэлектронных орбиталей:
Jastrow Function для спинов#
Дальше нас будут интересовать как раз модели Изинга и спины, а не частицы в пространстве или орбитали из вторичного квантования. Для спинов можем записать Jastrow function следующим образом:
где матрица \(W\) будет играть роль параметризации и отражать парные спиновые корреляции. Давайте посмотрим это на практике при помощи библиотеки NetKet [CCH+19] .
Моделировать будем простую модель Изинга для цепочки из 10 спинов (чтобы быстро считалось):
g = nk.graph.Hypercube(length=10, n_dim=1, pbc=True) hi = nk.hilbert.Spin(s=0.5, N=g.n_nodes) op = nk.operator.Ising(h=1.321, hilbert=hi, J=0.5, graph=g)
Поскольку модель относительно небольшая по числу частиц, то сразу можем получить точное решение методом Ланкзоша.
exact = nk.exact.lanczos_ed(op)[0]
Создадим модель на основе Jastrow и VMC:
sampler = nk.sampler.MetropolisLocal(hi) model = nk.models.Jastrow(dtype=complex) optimizer = nk.optimizer.Sgd(learning_rate=0.05) sr = nk.optimizer.SR(diag_shift=0.01) vmc = nk.driver.VMC(op, optimizer, sampler, model, n_samples=1008, preconditioner=sr)
/home/runner/work/qmlcourse/qmlcourse/.venv/lib/python3.8/site-packages/netket/utils/deprecation.py:126: FutureWarning: **DEPRECATION_WARNING:** The `dtype` argument to neural-network layers and models is deprecated throughout NetKet to maintain consistency with new releases of flax. Please use `param_dtype` instead. This warning will become an error in a future version of NetKet. warn_deprecation(_dep_msg) WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)
Изучение документации библиотеки NetKet оставляем вам самим, так как объяснение абстракций графа и гильбертова пространства, а также использование метода stochastic reconfiguration для вычисления градиентов выходит за рамки лекции. Документаци представлена на сайте NetKet.
logger = nk.logging.RuntimeLog() vmc.run(50, out=logger, show_progress=False)