=Paper=
{{Paper
|id=Vol-3041/584-591-paper-108
|storemode=property
|title=Параллельный Алгоритм Вычисления Функции Вигнера для Квантовой Системы с Полиномиальным Потенциалом (Parallel Algorithm for Calculating the Wigner Function for a Quantum System with a Polynomial Potential)
|pdfUrl=https://ceur-ws.org/Vol-3041/584-591-paper-108.pdf
|volume=Vol-3041
|authors=Evgeny Perepelkin,Boris Sadovnikov,Natalia Inozemtseva,Evgeny Burlakov,Rimma Polyakova,Pavel Sysoev,Marianna Sadovnikova
}}
==Параллельный Алгоритм Вычисления Функции Вигнера для Квантовой Системы с Полиномиальным Потенциалом (Parallel Algorithm for Calculating the Wigner Function for a Quantum System with a Polynomial Potential)==
https://ceur-ws.org/Vol-3041/584-591-paper-108.pdf
Proceedings of the 9th International Conference "Distributed Computing and Grid Technologies in Science and
Education" (GRID'2021), Dubna, Russia, July 5-9, 2021
PARALLEL ALGORITHM FOR CALCULATING
THE WIGNER FUNCTION FOR A QUANTUM SYSTEM
WITH A POLYNOMIAL POTENTIAL
E.E. Perepelkin1,2,4, B.I. Sadovnikov2, N.G. Inozemtseva3,4,
1, a
E.V. Burlakov2,4, R.V. Polyakova , P.N. Sysoev2, M.B. Sadovnikova2
1
Joint Institute for Nuclear Research, DUBNA, 141980 RUSSIA
2
Lomonosov Moscow State University, MOSCOW, 119991 RUSSIA
3
Dubna State University, DUBNA, 141980 RUSSIA
4
Moscow Technical University of Communications and Informatics, MOSCOW, 111024
RUSSIA
E-mail: a polykovarv@mail.ru
The article considers the construction of a parallel algorithm on the GPU computing architecture for
finding the Wigner function of a quantum system with a polynomial potential. A numerical-analytical
method for constructing the Wigner function, which is based on calculating the trace of the product of
the density matrix and the matrix of the Weyl operator, is described. The operators are represented in
the basis of a quantum harmonic oscillator, for which the Moyal equation transforms into the Liouville
equation. This approach enables to visually analyze the degree of anharmonicity of the system in terms
of the off-diagonal elements of the density matrix. The parallel implementation on the GPU massively
parallel architecture has reduced the calculation time by two orders of magnitude compared to the
single-threaded version on the x86 architecture.
Keywords: Wigner function, GPU computing, quantum system with a polynomial potential
Evgeny Perepelkin, Boris Sadovnikov, Natalia Inozemtseva, Evgeny Burlakov,
Rimma Polyakova, Pavel Sysoev, Marianna Sadovnikova
Copyright © 2021 for this paper by its authors.
Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
584
Proceedings of the 9th International Conference "Distributed Computing and Grid Technologies in Science and
Education" (GRID'2021), Dubna, Russia, July 5-9, 2021
ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ ВЫЧИСЛЕНИЯ
ФУНКЦИИ ВИГНЕРА ДЛЯ КВАНТОВОЙ СИСТЕМЫ
С ПОЛИНОМИАЛЬНЫМ ПОТЕНЦИАЛОМ
Е.Е. Перепёлкин1,2,4, Б.И. Садовников2, Н.Г. Иноземцева3,4,
Е.В. Бурлаков2,4, Р.В. Полякова1, a, П.Н. Сысоев2, М.Б. Садовникова2
1
Объединенный институт ядерных исследований (ОИЯИ), Дубна, Московская
область,
141980 Россия
2
Физический факультет МГУ имени М.В. Ломоносова, Москва, 119991 Россия
3
Университет «Дубна», Московская область,141980 Россия
4
Московский технический университет связи и информатики (МТУСИ), Москва,
111024 Россия
E-mail: a polykovarv@mail.ru
В работе рассмотрено построение параллельного алгоритма на вычислительной архитектуре
GPU нахождения функции Вигнера квантовой системы с полиномиальным потенциалом.
Описан численно-аналитический метод построения функции Вигнера, основанный на
вычислении следа произведения матрицы плотности и матрицы оператора Вейля.
Представление операторов выполнено в базисе квантового гармонического осциллятора, для
которого уравнение Моэля переходит в уравнение Лиувилля. Такой подход позволяет наглядно
анализировать степень ангармоничности системы по недиагональным элементам матрицы
плотности. Параллельная реализация на массивно-параллельной архитектуре GPU позволила
уменьшить время вычисления на два порядка в сравнении с однопоточной версией на x86
архитектуре.
Ключевые слова: функция Вигнера, вычисления на GPU, квантовая система с
полиномиальным потенциалом
Евгений Перепёлкин, Борис Садовников, Наталья Иноземцева, Евгений Бурлаков,
Римма Полякова, Павел Сысоев, Марианна Садовникова
Copyright © 2021 for this paper by its authors.
Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
585
Proceedings of the 9th International Conference "Distributed Computing and Grid Technologies in Science and
Education" (GRID'2021), Dubna, Russia, July 5-9, 2021
1. Введение
Решение современных задач квантовой информатики [1-3], квантовой связи и
криптографии [4, 5], задач обработки сигнала [6-8] содержит аппарат квантовой механики в
фазовом пространстве. Функция Вигнера [9] играет важную роль при нахождении средних по
фазовому пространству квантовых характеристик системы [10-20].
Эволюция функции Вигнера W r , p, t описывается уравнением Моэля
1 2 U , 2l 1 W . (i.1)
l 2l
W 1
p, r W rU , pW
t m l 1 2l 1! r p
где U — потенциал квантовой системы. В случае гармонического осциллятора
уравнение (i.1) имеет нулевую правую часть и совпадает с известным классическим
уравнением Лиувилля. Указанное упрощение уравнения (i.1) позволяет эффективно получать
выражения, описывающие эволюцию функции Вигнера. Для потенциалов степени выше второй
решение уравнения (i.1) затруднительно. В работе [21], используя уникальность потенциала
второй степени, были получены явные выражения для матричных элементов wn,k оператора
Вейля в базисе гармонического осциллятора ψ k .
W Tr ˆ Wˆ , (i.2)
s s ps
W= wn ,k x, p
det det 1
2
ψ x ψ x exp i
n k ds.
2 2
1 e x p P x i p , x i p ,
k 2
2 2
wn,k x, p
2 2
n ,k (i.3)
m
где k , n ck cn матрица плотности; . Полиномы Pn,k имеют вид:
min n , k
z1n s z2k s
Pn ,k z1 , z2 2 nk
n !k ! , (i.4)
s 0 2s s ! k s ! n s !
и удовлетворяют условию ортогональности с весом 2 x e x y :
2 2
x, y P x, y P
2 n1 , k1 n2 , k2 x, y dxdy N n 2n ,k k n n ,k k ,
1 2 1 2 1 2 1 2
(i.5)
min n1 , k1 min n2 , k2
k l s 1 !! n l s 1 !!
2
N n1 n2 ,k1 k2 n1 !n2 !k1 !k2 !
s 0
l 0 s ! n1 s ! k1 s ! l ! n2 l ! k2 l !
. где
s l , n , k even / odd
1, n, k чётные,
1, n, k нечётные,
n,k
0, n чётный, k нечётный,
0, n нечётный, k чётный.
При получении выражения (i.3) предполагалось, что волновая функция разлагается
по базису ψ k :
586
Proceedings of the 9th International Conference "Distributed Computing and Grid Technologies in Science and
Education" (GRID'2021), Dubna, Russia, July 5-9, 2021
x, t cn t ψn x . (i.6)
n 0
Отметим, что в этом случае выражение (i.2) может быть записано в виде свёртки:
W C T WC ,c (i.7)
C c1 , c2 ,... , C c1 , c2 ,... .
T T
Используя выражения (i.2)-(i.4), в работе [22] был описан метод построения функции
N
Вигнера для квантовой системы с полиномиальным потенциалом U N x an x n степени N .
n 1
Пусть волновая функция удовлетворяет уравнению Шрёдингера
2
xx U N E 0 и
2m
допускает разложение (i.6), тогда векторы коэффициентов C (i.7) являются собственными
векторами, а спектр энергий E= собственными значениями симметричной матрицы J n,k :
N
JC E=C , 0, J n,k I 0 al I nl,k ,
n n,k
(i.8)
l 2
l 1 l 2 n l 1 k l 1 n 0 k 0
I n ,k I n,k I n 1,k , I n 0,k n, k , I n1, k I n 1, k
l
I n,k 1 I n, k 1 ,
2 2 2 2 2
где m 2 2a2 ; m ;
1
n − спектр энергий гармонического осциллятора.
n
2
В данной работе описывается параллельный алгоритм, позволяющий эффективно
производить вычисление функции Вигнера на конечно-разностной сетке по приведенным выше
выражениям.
2. Параллельная реализация на графическом процессоре (GPU)
Для квантовой системы, характеризуемой волновыми функциями , функции
Вигнера допускают представление:
1 s
ps s
W x, p
2
x x exp i
ds, (1)
2 2
где — номер квантового состояния. Вычисления по формуле (1) требуют нахождения
порядка N p N x значений Wm , n на двумерной фазовой плоскости x, p , где N x и N p число
разбиений по координате и по импульсу соответственно ( m 0,..., N x ; n 0,..., N p ). При
численном интегрировании выражения (1) требуется «обрезание» пределов интегрирования.
Размер области интегрирования и количество разбиений N x и N p определяется характером
поведения волновых функций . Как правило, с увеличением номера состояния
увеличивается количество осцилляций волновой функции , что требует дополнительных
узлов разностной сетки при решении уравнения Шрёдингера и при численном интегрировании
выражения (1). Объём вычислений обусловлен не только нахождением сеточных функций Wm , n
587
Proceedings of the 9th International Conference "Distributed Computing and Grid Technologies in Science and
Education" (GRID'2021), Dubna, Russia, July 5-9, 2021
, но и необходимостью получения средних характеристик квантовой системы, требующих
повторного численного интегрирования с функцией Wm , n .
Ускорение вычислений сеточных функций Wm , n и средних характеристик квантовой
системы может быть достигнуто тремя способами.
Во-первых, использование явных выражений (i.3), (i.4) для матричных элементов
wn,k x, p позволяет избежать непосредственного численного интегрирования в выражении
(1). Функции wn,k x, p известны в явном виде и для фиксированной разностной сетки ( wn,k )
m,n
могут быть вычислены один раз и храниться в памяти.
Во-вторых, благодаря возможности разложения волновой функции по базису
гармонического осциллятора ψ j (i.6):
cj ψ j , (2)
j 0
1
m x
m 4 2 m
2
1
ψ j x e H j x , j 0 , (3)
2 j!
0
j
можно с одной стороны, вместо количества разбиений сеточной функции контролировать
количество членов N C, ряда (2), а с другой стороны, быстрый экспоненциальный спад (3)
гарантирует небольшое значение N C, для приемлемой точности аппроксимации волновой
функции конечным разложением по функциям ψ j (2). Согласно (i.7), (i.8) коэффициенты
разложения c j образуют собственные вектора C матрицы J , вид которой известен точно
(i.8) и может быть вычислен один раз до начала основной процедуры вычислений. Отметим,
что матрица J является симметричной, что уменьшает размер памяти для ее хранения и
позволяет воспользоваться специальными численными методами при решении задачи на
собственные значения.
В-третьих, задача на собственные значения (поиск C и E=), вычисление значений
wn m,k,n и свертки (i.7) Wm ,n , а также нахождение квантовых средних по функции Wm ,n являются
процедурами, допускающими эффективное численное распараллеливание. Подходящей
вычислительной архитектурой для таких задач может быть массивно-параллельная архитектура
графических процессоров GPU с использованием среды CUDA [23, 24]. Задача на собственные
значения, операции по перемножению матриц имеют реализацию в бесплатных пакетах
библиотек компании NVIDIA: cuBLAS, cuSPARS, cuSolver.
В качестве примера рассмотрим алгоритм параллельного вычисления матричных
m,n
элементов wn,k . Матричные элементы wn,k обладают свойством wn,k wk ,n , поэтому в
памяти достаточно хранить только диагональ и верхнюю треугольную часть матрицы W .
Размер матрицы W определяется числом N C, , которое удовлетворяет условию NC, N x или
NC, N p . Количество разбиений N x и N p можно считать величинами одного порядка:
588
Proceedings of the 9th International Conference "Distributed Computing and Grid Technologies in Science and
Education" (GRID'2021), Dubna, Russia, July 5-9, 2021
n,m n,m
w1,1 w1,2 ... w1, nN,m
C,
n,m n ,m n ,m
W = 0 w2,2 .. w3, N . (4)
C,
... 0 ...
Заметим, что при желании можно выбрать максимальное значение квантового
состояния L max (для которых будут производиться расчёты) и задать размерность
матрицы (4) NC , L NC , L .
Так как величина N x N p определяет максимальный объём вычислений, то логично
произвести распараллеливание именно по индексам n, m , задав двумерную сетку нитей и
блоков nx, np . Заголовок функции-ядра (kernel-function) имеет вид:
__global__ void W_Matrix (int n, int k, double *WM...)
{int nx = threadIdx.x + blockIdx.x * blockDim.x;
int np = threadIdx.y + blockIdx.y * blockDim.y;
// (x, p) parallelism ...
}
Каждая вычислительная нить (thread) находит свой матричный элемент в
соответствующей ее номеру точке фазового пространства (см. рис. 1). Функция-ядро
2
«W_Matrix» запускается порядка 12 NC , L раз c «host», чтобы произвести вычисления всех
элементов матрицы (4). Результаты вычислений можно хранить в глобальной памяти GPU
(global memory) с целью их переиспользования другими функциями-ядрами.
Рис. 1. Распределение вычислительной нагрузки по потокам
3. Заключение
Время работы параллельной реализации алгоритма (GPU GTX1050) вычисления
функции Вигнера на два порядка меньше времени выполнения однопоточной версии (CPU
класса Core i7).
4. Благодарность
Работа выполнена при поддержке гранта РФФИ No. 18-29-10014 и Междисциплинарной
научно-образовательной школы МГУ им. М.В. Ломоносова «Фотонные и квантовые
технологии. Цифровая медицина»
589
Proceedings of the 9th International Conference "Distributed Computing and Grid Technologies in Science and
Education" (GRID'2021), Dubna, Russia, July 5-9, 2021
Литература
[1] R. P. Rundle, Todd Tilma, J. H. Samson, V. M. Dwyer, R. F. Bishop, and M. J. Everitt, General
approach to quantum mechanics as a statistical theory, Phys. Rev. A 99, 012115 – Published 16
January 2019
[2] Ievgen I.Arkhipov, Artur Barasiński, Jiří Svozilík, Negativity volume of the generalized Wigner
function as an entanglement witness for hybrid bipartite states, Sci Rep 8, 16955 (2018)
[3] Andersen, U., Neergaard-Nielsen, J., van Loock, P. et al. Hybrid discrete- and continuous-
variable quantum information. Nature Phys 11, 713–719 (2015).
[4] Alberto Casado, Santiago Guerra, José Plácido, From Stochastic Optics to the Wigner
Formalism: The Role of the Vacuum Field in Optical Quantum Communication Experiments, Atoms
2019, 7, 76
[5] Casado, A.; Guerra, S.; Plácido, J. Wigner representation for experiments on quantum
cryptography using two-photon polarization entanglement produced in parametric down-conversion. J.
Phys. B At. Mol. Opt. Phys. 2008, 41, 045501.
[6] Cohen L. Time-Frequency Analysis. Prentice Hall, Englewood Cliffs, 1995
[7] Zayed, A. A New Perspective on the Two-Dimensional Fractional Fourier Transform and Its
Relationship with the Wigner Distribution . J Fourier Anal Appl 25, 460–487 (2019)
[8] Claasen, T.A.C.M., Mecklenbräuker,W.F.G.: TheWigner distribution—a tool for time-
frequency signal analysis. II: discrete-time signals, part 2. Philips J. Res. 35, 276–300 (1980)
[9] E.P. Wigner, On the quantum correction for thermodynamic equilibrium, Phys. Rev. 40 (June
1932) 749—759
[10] Meystre P., Walls D.F. Nonclassical Effects in Quantum Optics. American Institute of Physics,
New York, 1991
[11] Kimble H.J., Walls D.F. // J. Opt. Soc. Am. B. 1987. V. 4A0)
[12] Knight P., Loudon R. // J. Mod. Opt. 1984. V. 34F)
[13] Giacobino E., Fabre С // Appl. Phys. B. 1992. V. 55C)
[14] Wu L.A., Xiao М., Kimble H.J. Squeezed states of light from an optical parametric oscillator //
J. Opt. Soc. Am. B. 1987. V. 4. P. 1465-1475.
[15] Schleich W.P., Wheeler J.A. Oscillations in the Photon Distribution of Squeezed States and
Interference in Phase Space // Nature (London) 326, 574-577 A987a)
[16] Smithey D.T., Beck M., Raymer M.G., Faridani A. Measurement of the Wigner distribution and
the density matrix of a light mode using optical homodyne tomography: application to squeezed states
and the vacuum // Phys. Rev. Lett. 1993. V. 70. P. 1244-1247.
[17] Breitenbach G., Schiller S., Mlynek J. Measurement of the quantum states of squeezed light //
Nature. 1997. V. 387. P. 471-475.
[18] Royer A. Measurement of the Wigner function // Phys. Rev. Lett. 1985. V. 55. P. 2745-2748.
[19] Vogel K., Risken H. Determination of quasiprobability distributions in terms of probability
distributions for the rotated quadrature phase // Phys. Rev. A. 1989. V. 40. P. 2847-2849.
[20] Radon J. fiber die Bestimmung von Funktionen durch ihre Integralwerte langs gewisser
Mannigfaltigkeiten // Ber. Verh. Sachs. Akad. Wiss. Leipzig, Math.-Nat. Kl. 1917. V.69. P. 262-277.
590
Proceedings of the 9th International Conference "Distributed Computing and Grid Technologies in Science and
Education" (GRID'2021), Dubna, Russia, July 5-9, 2021
[21] Perepelkin E.E., Sadovnikov B.I., Inozemtseva N.G., Burlakov E.V., Explicit form for the
kernel operator matrix elements in eigenfunction basis of harmonic oscillator, Journal of Statistical
Mechanics: Theory and Experiment, (2020) № 023109
[22] Perepelkin E.E., Sadovnikov B.I., Inozemtseva N.G., Burlakov E.V., Wigner function of a
quantum system with polynomial potential, Journal of Statistical Mechanics: Theory and Experiment,
(2020) № 053105
[23] Н.Г. Иноземцева, Е.Е. Перепёлкин, Б.И. Садовников, Оптимизация алгоритмов задач
математической физики для графических процессоров, Физический факультет МГУ имени
М.В. Ломоносова Москва, 2012, ISBN 978-5-8279-0107-5, 256 с.
[24] Е.Е. Перепёлкин, Б.И. Садовников, Н.Г. Иноземцева, Вычисления на графических
процессорах (GPU) в задачах математической и теоретической физики, Серия «Классический
учебник МГУ», URSS Москва, 2019, ISBN 978-5-9710-6490-9, 240 с.
591