Однокубитное преобразование вектора-состояния

Однокубитное преобразование вектора-состояния
Последовательный алгоритм
Последовательная сложность	[math]3 \cdot 2^n[/math]
Объём входных данных	[math]2^n+4[/math]
Объём выходных данных	[math]2^n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]2[/math]
Ширина ярусно-параллельной формы	[math]2^{n+1}[/math]

Основные авторы описания: А.Ю.Чернявский, Вад.В.Воеводин (раздел 2.2)

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

Алгоритм производит моделирования действия однокубитного квантового вентиля на вектор-состояние. Данный алгоритм обычно является подпрограммой и многократно применяется к различным кубитам одного состояния (например при моделировании квантовых алгоритмов или анализе квантовой запутанности). Особенностью алгоритма, как и большинства алгоритмов квантовой инфоорматики, является экспоненциальный рост объема данных в зависимости от основного параметра - числа кубитов, что приводит к необходимости суперкомпьютерной реализации для решения важных практических задач.

1.2 Математическое описание алгоритма

Исходные данные:

Целочисленные параметры [math]n - [/math] число кубитов (необязательно) и [math]k -[/math] номер кубита, над которым производится преобразование.

Комплексная матрица [math]U = \begin{pmatrix} u_{00} & u_{01}\\ u_{10} & u_{11} \end{pmatrix}[/math] однокубитного преобразования размера [math]2 \times 2.[/math]

Комплексный вектор [math]v[/math] размерности [math]2^n,[/math] задающей начальное состояние многокубитной системы.

Вычисляемые данные: комплексный вектор [math]w[/math] размерности [math]2^n,[/math] соответствующий состоянию после преобразования.

Формулы метода:

Состояние после действия преобразования [math]U[/math] на [math]k-[/math]й кубит имеет вид [math]v_{out} = I_{2^{k-1}}\otimes U \otimes I_{2^{n-k}},[/math] где [math]I_{j} - [/math] единичная матрица размерности [math]j,[/math] а [math]\otimes - [/math] тензорное произведение (произведение Кронекера).

Однако, элементы итогового вектора можно записать и в прямом виде, что более удобно для вычислений:

[math] w_{i_1i_2\ldots i_k \ldots i_n} = \sum\limits_{j_k=0}^1 u_{i_k j} v_{i_1i_2\ldots j_k \ldots i_n} = u_{i_k 0} v_{i_1i_2\ldots 0_k \ldots i_n} + u_{i_k 1} v_{i_1i_2\ldots 1_k \ldots i_n} [/math]

Индекс-кортеж [math]i_1i_2\ldots i_n[/math] представляет собой двоичную запись индекса элемента в массиве.

1.3 Вычислительное ядро алгоритма

Вычислительное ядро алгоритма представляет собой независимое вычисление всех [math]2^n[/math] элементов вектора [math]w.[/math] Вычисление каждого элемента требует две операции умножения и одну операцию сложения. Кроме того необходимо вычислять индексы типа [math]i_1i_2\ldots 0_k \ldots i_n,[/math] а также значение бита [math]i_k,[/math] что требует побитовых операций.

1.4 Макроструктура алгоритма

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

1.5 Схема реализации последовательного алгоритма

Для индекса [math]i[/math] от [math]0[/math] до [math]2^n-1[/math]

1. Вычислить элемент [math]i_k[/math] двоичного представления индекса [math]i.[/math]
2. Вычислить индексы [math]j[/math] имеющие двоичные представления [math]i_1i_2\ldots \overline{i_k} \ldots i_n,[/math] где крышка означает обращение бита.
3. Вычислить [math]w_i = u_{i_k i_k}\cdot v_{i} + u_{i_k \overline{i_k}}\cdot v_j.[/math]

1.6 Последовательная сложность алгоритма

Алгоритм требует:

1. [math]2^{n+1}[/math] операций умножения комплексных чисел;
2. [math]2^n[/math] операций сложения комплексных чисел;
3. [math]2^n[/math] операций получения значения [math]k[/math]-го бита числа;
4. [math]2^n[/math] операций изменения значения [math]k[/math]-го бита числа.

Отметим, что данный алгоритм обычно применяется много раз подряд, в связи с чем вычисления, связанные с побитовыми операциями (3-4), могут единожды проводиться в начале алгоритма. Кроме того, от них можно избавиться, пользуясь сложением и логическим умножением с числом [math]2^k,[/math] которое сохраняется для всего алгоритма.

1.7 Информационный граф

Представим рисунки графов алгоритма для случая [math]n=3, k=1[/math] (рис.1) и [math]k=2[/math] (рис.2). На графах не представлены матрицы преобразования [math]U,[/math] в связи с тем, что их размер при больших [math]n[/math] много меньше, нежели размеры входного и выходного векторов. На Рис.3 изображена основная операция, представляемая на Рис.1 и Рис.2 оранжевым.

Отображение графа со входными и выходными данных, а также "свёрнутой" тройной операцией удобно для понимания локальности обращений к памяти. Отметим, что структура графа (а именно обращение к входным данным) сильно зависит от параметра [math]k.[/math]

Рисунок 1. Граф алгоритма для [math]n=3, k=1[/math] без отображения матрицы преобразования [math]U.[/math]

Рисунок 2. Граф алгоритма для [math]n=3, k=2[/math] без отображения матрицы преобразования [math]U.[/math]

Рисунок 3. Граф основной операции.

1.8 Ресурс параллелизма алгоритма

Как видно из информационного графа, прямой алгоритм моделирования однокубитного преобразования обладает высочайшей степенью параллелизма. Все операции вычисления элементов нового вектора-состояния могут быть произведены параллельно. Для вычисления одного элемента необходимо выполнить две операции умножения и одну операцию сложения, операции умножения, в свою очередь, также могут быть выполнены параллельно. Таким образом, необходимо выполнение двух ярусов: одного, состоящего из [math]2^{n+1}[/math] умножений и другого, состоящего из [math]2^n[/math] сложений.

1.9 Входные и выходные данные алгоритма

Входные данные:

• Комплексный вектор состояния [math]u[/math] длины [math]2^n.[/math] Обычно нормирован на единицу.
• Унитарная матрица [math]U[/math] порядка [math]2[/math].
• Номер кубита [math]q[/math].

Выходные данные:

• Комплексный вектор состояния [math]w[/math] длины [math]2^n[/math].

Объем входных данных: [math]2^n+4[/math] комплексных чисел и [math]1[/math] целочисленный параметр.

Объем выходных данных: [math]2^n[/math] комплексных чисел.

1.10 Свойства алгоритма

Соотношение последовательной и параллельной сложности является экспоненциальным (эксмпонента переходит в константу).

Вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных – константа.

2 Программная реализация алгоритма

2.1 Особенности реализации последовательного алгоритма

На языке C функцию однокубитного преобразования можно записать следующим образом:

void OneQubitEvolution(complexd *in, complexd *out, complexd U[2][2], int nqubits, int q)
{
	//n - число кубитов
        //q - номер кубита для преобразования

	int shift = nqubits-q;
	//Все биты нулевые, кроме соответствующего позиции преобразуемого кубита  
	int pow2q=1<<(shift);

	int N=1<<nqubits;
	for	(int i=0; i<N; i++)
	{
		//Обнуления меняющегося бита
		int i0 = i & ~pow2q;
		
		//Установка меняющегося бита
		int i1 = i | pow2q;
		
		//Получение значения меняющегося бита
		int iq = (i & pow2q) >> shift;

		out[i] = U[iq][0] * in[i0] + U[iq][1] *in[i1];
	}
}

Отметим, что существенная часть вычислений и логики кода приходится на битовые операции, однако, этого можно избежать: однокубитное преобразование в большинстве случаев является лишь подпрограммой и применяется к разным кубитам большое число раз. В свою очередь, вычисляемые при помощи битовых операций индексы i0, i1 и iq зависят лишь от параметров количества кубитов n и номера кубита q. Число кубитов обычно фиксировано, соответственно, можно вычислить эти индексы заранее для всех q от 1 до n. Для хранения потребуется лишь массив целочисленных переменных линейного размера 3n в то время, как обрабатываемые данные имеют экспоненциальный размер. Очевидно, что такая оптимизация критически необходима при реализации алгоритма на вычислительных устройствах или языках программирования с отсутствием быстрых битовых операций (примером может служить среда Matlab).

2.2 Локальность данных и вычислений

К сожалению, в потивовес идельной возможности параллелизации алгоритма, он обладает очень плохой локальностью.

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

Однако, преобразование одного кубита в прикладных задачах является лишь подпрограммой и применяется многократно с различными параметрами q. Как видно из математического описания и Рис. 1-2, это полностью исключает возможность добиться локальности обращений к данным.

2.2.1 Локальность реализации алгоритма

2.2.1.1 Структура обращений в память и качественная оценка локальности

Рисунок 1. Однокубитное преобразование вектора-состояния. Общий профиль обращений в память

На рис. 1 представлен профиль обращений в память для вычисления однокубитного преобразования вектора-состояния. Данный профиль состоит из обращений к трем массивам, фрагменты для отдельных массивов выделены на рис. 1 зеленым цветом. Из общего профиля можно увидеть, что обращения редко используются повторно, по крайней мере в случае фрагментов 2 и 3. При этом обращения к близко расположенным друг к другу данным выполняются рядом. Рассмотрим выделенные фрагменты поближе.

Отдельно фрагмент 1 представлен на рис. 2. Видно, что данный массив состоит всего из 4-х элементов, к которым постоянно выполняются обращения. Такой фрагмент обладает очень высокой локальностью, поскольку постоянно переиспользуются запрошенные ранее данные.

Далее рассмотрим фрагмент 2 (рис. 3). Здесь все еще проще – выполняется обычный последовательный перебор всех элементов массива. Такой фрагмент обладает высокой пространственной локальностью, однако очень низкой временной (данные просто не используются повторно).

Рисунок 2. Фрагмент 1 (профиль обращений к первому массиву)

Рисунок 3. Фрагмент 2 (профиль обращений ко второму массиву)

Наиболее интересным представляется фрагмент 3. Его небольшой фрагмент, выделенный на рис. 1 желтым, представлен на рис. 4. Однако при ближайшем рассмотрении оказывается, что данный фрагмент тоже просто устроен, хотя и немного сложнее предыдущих.

В данном случае также виден в центре последовательный перебор всех элементов массива, параллельно с которым выполняются обращения либо к элементам с большим или меньшим виртуальным адресом. Отметим, однако, что эта разница между виртуальными адресами, судя по всему, больше 64 байт (длины строки), что может служить причиной возникновения большого числа кэш-промахов.

Рисунок 4. Небольшая часть фрагмента 3, выделенная на рис. 1 желтым

В общем можно сказать, что общий профиль обращений в память обладает достаточно высокой пространственной локальностью, поскольку большинство обращений образуют последовательные переборы элементов массивов, однако временная локальность низка – данные практически не используются повторно.

2.2.1.2 Количественная оценка локальности

Основной фрагмент реализации, на основе которого были получены количественные оценки, приведен здесь (функция Kernel). Условия запуска описаны здесь.

Первая оценка выполняется на основе характеристики daps, которая оценивает число выполненных обращений (чтений и записей) в память в секунду. Данная характеристика является аналогом оценки flops применительно к работе с памятью и является в большей степени оценкой производительности взаимодействия с памятью, чем оценкой локальности. Однако она служит хорошим источником информации, в том числе для сравнения с результатами по следующей характеристике cvg.

На рисунке 5 приведены значения daps для реализаций распространенных алгоритмов, отсортированные по возрастанию (чем больше daps, тем в общем случае выше производительность). Можно увидеть, что производительность работы с памятью для этой программы высока – значение daps примерно на уровне теста Linpack. Видимо, низкая временная локальность в данном случае компенсируется высокой пространственной локальностью.

Рисунок 5. Сравнение значений оценки daps

Вторая характеристика – cvg – предназначена для получения более машинно-независимой оценки локальности. Она определяет, насколько часто в программе необходимо подтягивать данные в кэш-память. Соответственно, чем меньше значение cvg, тем реже это нужно делать, тем лучше локальность.

На рисунке 6 приведены значения cvg для того же набора реализаций, отсортированные по убыванию (чем меньше cvg, тем в общем случае выше локальность). Можно увидеть, что, в отличие от daps, cvg оценивает локальность данной программы как достаточно низкую. В частности, значение cvg для Linpack заметно меньше, в то время как значения daps практически совпадали.

Одна из возможных причин этого – влияние арифметических операций. Может получиться, что данные из памяти не будут запрашиваться, пока арифметические операции не будут выполнены; это приводит к простою подсистемы памяти. Соответственно, если в одной программе таких операций нет, а в другой - есть, то daps в первом случае будет выше. При этом cvg не поменяется, поскольку эта оценка не зависит от времени выполнения.

В данном случае арифметических операций практически нет (в отличие от некоторых других программ), поэтому daps может показывать более высокие результаты, в то время как cvg показывает достаточно низкую оценку.

Рисунок 6. Сравнение значений оценки cvg

2.3 Возможные способы и особенности параллельной реализации алгоритма

Основной способ параллельной реализации очевиден - необходимо распараллеливание основного цикла (параллельное вычисление различных компонент выходного вектора-состояния) и, желательно, операций умножения. На машинах с общей памятью такой вариант распараллеливания приводит к ускорению, близкому к максимально-возможному. Однако, данный способ сталкивается с проблемами локальности данных. Анализируя математическая описание и информационные графы алгоритма легко видеть, что при использовании большого числа узлов с собственной памятью количество необходимых пересылок между узлами становится сопоставимым с количеством вычислений, что приводит к существенной потере эффективности. Возможны разные пути частичного решения этой проблемы, например, кэширование или использование парадигмы программирования SHMEM, однако, столь сильная нелокальность использования данных всё равно не позволяет получить хорошей эффективности.

2.4 Масштабируемость алгоритма и его реализации

2.4.1 Масштабируемость алгоритма

В связи с констатной высотой ярусно-параллельной формы алгоритм имеет отличную масштабируемость на машинах с общей памятью. При реализации на машинах с разделяемой памятью быстро возникают проблемы, связанные с локальностью.

Набор и границы значений изменяемых параметров запуска реализации алгоритма:

• число процессоров [4 : 512] со значениями степени двойки;
• число кубитов [16;32] с шагом 2.

В результате проведённых экспериментов был получен следующий диапазон эффективности реализации алгоритма:

• минимальная эффективность реализации 0.000003%;
• максимальная эффективность реализации 0.008%.

На следующих рисунках приведены графики производительности и эффективности выбранной реализации алгоритма в зависимости от изменяемых параметров запуска.

Рисунок 8. Параллельная реализация алгоритма. Изменение производительности в зависимости от числа процессоров и числа кубитов.

Рисунок 9. Параллельная реализация алгоритма. Изменение эффективности в зависимости от числа процессоров и числа кубитов.

2.6 Выводы для классов архитектур

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