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

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

Основные авторы описания: А.Ю.Чернявский

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] комплексных чисел [math]1.[/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.3 Возможные способы и особенности параллельной реализации алгоритма

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

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

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

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