Single-qubit transform of a state vector

Single-qubit transform of a state vector
Sequential algorithm
Serial complexity	[math]3 \cdot 2^n[/math]
Input data	[math]2^n+4[/math]
Output data	[math]2^n[/math]
Parallel algorithm
Parallel form height	[math]2[/math]
Parallel form width	[math]2^{n+1}[/math]

Primary authors of this description: A.Yu.Chernyavskiy, Vad.V.Voevodin (Section 2.2)

1 Properties and structure of the algorithm

1.1 General description of the algorithm

The algorithm simulates an action of a single-qubit quantum gate on a state vector.

This algorithm is often a subroutine and is used iteratively to the different qubits of a state. For example, it's used for the quantum algorithms simulation and the analysis of quantum entanglement. As for the most of the classical algorithms in quantum information science, the described algorithm shows the exponential data growth with respect to the number of qubits. This fact leads the necessity of a supercomputer usage for the practically important tasks.

1.2 Mathematical description of the algorithm

Input data:

Integer parameters [math]n[/math] (the number of qubits) and [math]k[/math] (the number of the affected qubit).

A complex [math]2 \times 2[/math] matrix [math]U = \begin{pmatrix} u_{00} & u_{01}\\ u_{10} & u_{11} \end{pmatrix}[/math] of a single-qubit transform.

A complex [math]2^n[/math]-dimensional vector [math]v,[/math] which denotes the initial state of a quantum system.

Output data: the complex [math]2^n[/math]-dimensional vector [math]w,[/math] which corresponds to the quantum state after the transformation.

The single-qubit transform can be represented as follows:

The output state after the action of a transform [math]U[/math] on the [math]k-[/math]th qubit is [math]v_{out} = I_{2^{k-1}}\otimes U \otimes I_{2^{n-k}},[/math] where [math]I_{j}[/math] is the [math]j-[/math]dimensional identity matrix, and [math]\otimes[/math] denotes the tensor (Kronecker) product.

However, the elements of the output vector can be represented in more computationally useful direct form:

[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]

A tuple-index [math]i_1i_2\ldots i_n[/math] denotes the binary form of an array index.

1.3 The computational kernel of the algorithm

The computational kernel of the single-qubit transform can be composed of [math]2^n[/math] calculations of the elements of the output vector [math]w.[/math] The computation of each element consists of two complex multiplication and one complex addition opearations. Moreover we need to compute an indexes [math]i_1i_2\ldots 0_k \ldots i_n,[/math] and the [math]i_k[/math] bits, these computations needs the bitwise operations.

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

As noted in the description of the computational kernel, the main part of the algorithm consists of the independent computations of the output vector elements.

1.5 Implementation scheme of the serial algorithm

For [math]i[/math] from [math]0[/math] to [math]2^n-1[/math]

1. Calculate the element [math]i_k[/math] of the binary representation of the index [math]i.[/math]
2. Calculate indexes [math]j[/math] with binary representation [math]i_1i_2\ldots \overline{i_k} \ldots i_n,[/math] where the overline denotes the bit-flip.
3. Calculate [math]w_i = u_{i_k i_k}\cdot v_{i} + u_{i_k \overline{i_k}}\cdot v_j.[/math]

1.6 Serial complexity of the algorithm

The following number of operations is needed for the single-qubit transform:

1. [math]2^{n+1}[/math] complex multiplications;
2. [math]2^n[/math] complex additions;
3. [math]2^n[/math] [math]k[/math]-th bit reads;
4. [math]2^n[/math] [math]k[/math]-th bit writes.

Note, that the algorithm are often used iteratively, that's the the bit operations (3-4) can be performed once in the beginning of computations. Moreover, they can be avoided by the usage of additions and bitwise multiplications with the [math]2^k,[/math] that can be also calculated once for each [math]k.[/math]

1.7 Information graph

The Figs. 1 and 2 shows the information graph for the paramters [math]n=3, k=1[/math] and [math]n=3, k=2.[/math] The graphs are presented without a transformation matrix [math]U,[/math] because its small size in comparison with a input and output data. Fig.3 shows the main operation, which is denoted by orange on Figs.1-2.

The representation with an input and data and the "folded" triple operation is useful for the understanding of the memory locality. Note, the the structure of an information graph strongly depends on the parameter [math]k.[/math]

Figure 1. The information graph for [math]n=3, k=1[/math] without transformation matrix [math]U.[/math]

Figure 2. The information graph for [math]n=3, k=2[/math] without transformation matrix [math]U.[/math]

Figure 3. The information graph of the main operation.

1.8 Parallelization resource of the algorithm

As we can see on the information graph, the direct algorithm of one-qubit transform simulation has the very high parallelization resource, because all computations of different output vector elements can be made independently. Two multiplication operations needed for the computation of each output vector element also can be done in parallel.

So, there are only two levels of computations:

1. [math]2^{n+1}[/math] multiplications.
2. [math]2^n[/math] additions.

1.9 Input and output data of the algorithm

Input data:

• A [math]2^n[/math]-dimensional complex state vector [math]u.[/math] In most cases thу vector is normalized to 1.
• A [math]2[/math]-dimensional unitary matrix [math]U[/math]>.
• The number (index) of the affected qubit [math]q[/math].

Output data:

• The [math]2^n[/math]-dimensional complex output state vector [math]w.[/math].

Amount of input data: [math]2^n+4[/math] complex numbers and [math]1[/math] integer parameter.

Amount of output data: [math]2^n[/math] complex numbers.

1.10 Properties of the algorithm

The ratio of the serial complexity to the parallel complexity is exponential (an exponent transforms to a constant).

The computational power of the algorithm considered as the ratio of the number of operations to the amount of input and output data is constant.

2 Software implementation of the algorithm

2.1 Implementation peculiarities of the serial algorithm

The function of the single-qubit transform may be represented in C language as:

void OneQubitEvolution(complexd *in, complexd *out, complexd U[2][2], int nqubits, int q)
{
	//n - number of qubits
        //q - index of the affected qubit

	int shift = nqubits-q;
	//All bits are zeros, except the bit correspondent to the affected qubit.  
	int pow2q=1<<(shift);

	int N=1<<nqubits;
	for	(int i=0; i<N; i++)
	{
		//Making the changing bit zero
		int i0 = i & ~pow2q;
		
		//Setting the changing bit value
		int i1 = i | pow2q;
		
		//Getting the changing bit value
		int iq = (i & pow2q) >> shift;

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

Note, that the huge part of the computation and code logic consists of the bit operations, but this can be avoided. In most cases, single-qubit transform is just as a subroutine and is used sequentially many times. Moreover, indexes i0, i1, and iq depends only on the parameters n and q. The n is fixed, and we can calculate these indexes in the beginning of computations for every parameter q, and we need only 3n integers to store the data, when the input and output data growths exponentially. It's obvious that such optimization is critically needed for the implementations of the algorithm for the hardware and software with slow bit operations (for example, 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.2 Масштабируемость реализации алгоритма

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

Проведём исследование масштабируемости параллельной реализации алгоритма согласно методике. Исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета. Набор и границы значений изменяемых параметров запуска реализации алгоритма:

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

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

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

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

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

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

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

Исследованная параллельная реализация алгоритма

2.5 Динамические характеристики и эффективность реализации алгоритма

Для проведения экспериментов использовалась реализация алгоритма однокубитного преобразования, в реализации, доступной здесь. Все результаты получены на суперкомпьютере "ГрафИТ!". Использовались процессоры Intel Xeon X5570 с пиковой производительностью в 94 Гфлопс, а также компилятор Gnu 4.4.7. На рисунках показана эффективность реализации алгоритма встречной прогонки.

Рисунок 9. График загрузки CPU при выполнении алгоритма однокубитного преобразования

На графике загрузки процессора видно, что почти все время работы программы уровень загрузки составляет около 50% в среднем. Это указывает на высокую нагрузку выполнения алгоритма при котором на 8-ми ядерном процессоре работают 8 процессов, без использования технологии HyperTreading. Это означает, что использование всех физических ядер близко к 100%

Рисунок 10. График операций с плавающей точкой в секунду при выполнении алгоритма однокубитного преобразования

На Рисунке 11 показан график количества операций с плавающей точкой в секунду. На графике видна общая достаточно высокая производительность вычислений, достигающая в пике 300 млн операций в секунду и около 70 млн операций в секунду в среднем по всем узлам. Это может указывать на некоторый дисбаланс в вычислениях, при общем сравнительно высоком уровне производительности.

Рисунок 11. График кэш-промахов L1 в секунду при работе алгоритма однокубитного преобразования

На графике кэш-промахов первого уровня видно, что число промахов достаточно высокое и находится на уровне 4 млн/сек для отдельных процессов и 1 млн/сек в среднем по всем процессам. Это указывает на достаточно высокую интенсивность вычислений и хорошую локальность данных, учитывая высокие показатели операций с плавающей точкой и на дисбаланс вычислений между процессами.

Рисунок 12. График кэш-промахов L3 в секунду при работе алгоритма однокубитного преобразования

На графике кэш-промахов третьего уровня видно, что число промахов все еще достаточно высокое для параллельного приложения и находится на уровне 1 млн/сек в пике и около 250 тыс/сек в среднем по всем узлам. Отношение промахов L1/L3 около 4, что является средним показателем по таким задачам. Это также указывает на достаточно плохую локальность вычислений.

Рисунок 13. График количества чтений из оперативной памяти при работе алгоритма однокубитного преобразования

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

Рисунок 14. График количества записей в оперативную память при работе алгоритма однокубитного преобразования

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

Рисунок 15. График числа процессов, ожидающих вхождения в стадию счета (Loadavg), при работе алгоритма однокубитного преобразования

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

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

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