Уровень алгоритма

Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы, вещественный точечный вариант

Материал из Алговики
Перейти к навигации Перейти к поиску


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

Основные авторы описания: А.В.Фролов

Содержание

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

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

Метод Хаусхолдера (в советской математической литературе чаще называется методом отражений) используется для разложения матриц в виде [math]A=QR[/math] ([math]Q[/math] - унитарная, [math]R[/math] — правая треугольная матрица)[1]. При этом матрица [math]Q[/math] хранится и используется не в своём явном виде, а в виде произведения матриц отражения[2]. Каждая из матриц отражения может быть определена одним вектором. Это позволяет в классическом исполнении метода отражений хранить результаты разложения на месте матрицы A с использованием минимального одномерного дополнительного массива.

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

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

В методе Хаусхолдера для выполнения [math]QR[/math]-разложения матрицы используются умножения слева её текущих модификаций на матрицы Хаусхолдера (отражений).

Матрица отражений (Хаусхолдера) - матрица вида [math]U=E-2ww^*[/math], где [math]w[/math] - вектор, удовлетворяющий равенству [math]w^{*}w=1[/math]. Является одновременно унитарной ([math]U^{*}U=E[/math]) и эрмитовой ([math]U^{*}=U[/math]), поэтому обратна самой себе ([math]U^{-1}=U[/math]).

На [math]i[/math]-м шаге метода с помощью преобразования отражения "убираются" ненулевые поддиагональные элементы в [math]i[/math]-м столбце. Таким образом, после [math]n-1[/math] шагов преобразований получается матрица [math]R[/math] из [math]QR[/math]-разложения.

На каждом из шагов метода матрицу отражений обычно представляют не в стандартном виде, а в виде [math]U=E-\frac{1}{\gamma}vv^*[/math], где [math]v[/math] находится через координаты текущего [math]i[/math]-го столбца так:

[math]s[/math] - вектор размерности [math]n+1-i[/math], составленный из элементов [math]i[/math]-го столбца, начиная с [math]i[/math]-го.

Если [math](s,s)=0[/math], то [math]v=e_{i}[/math], [math]\gamma = \frac{1}{2}[/math].

В остальных случаях по алгоритму вычисляется [math]u = \frac{1}{\sqrt{(s,s)}}s[/math], и далее [math]v_{j}=0[/math] при [math]j\lt i[/math], [math]v_{j}=u_{j-i+1}[/math] при [math]j\gt i[/math], а [math]v_{i}=1[/math], если [math]u_{1}=0[/math] и [math]v_{i}=\frac{u_{1}}{|u_{1}|}(1+|u_{1}|)[/math] для остальных значений. При этом [math]\gamma =1+|u_{1}|=|v_{i}|[/math].

После вычисления вектора [math]v[/math] подстолбцы справа от ведущего модифицируются по формулам [math]x'=x-\frac{(x,v)}{\gamma}v[/math].

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

Основную часть алгоритма составляют вычисления на каждом шагу скалярных произведений [math](s,s)[/math] и [math](x,v)[/math] для всех подстолбцов [math]x[/math]справа от текущего, а также векторные операции [math]x'=x-\frac{(x,v)}{\gamma}v[/math]. Это используется при программировании метода во многих библиотеках для его конструирования из стандартных подпрограмм (например, из BLAS).

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

Рисунок 1. Граф шага (обнуление [math]i[/math]го столбца) алгоритма. Квадратики - входные данные шага (берутся с предыдущего или из входных данных), кружки - операции. Зелёным выделены операции типа a+bb, салатовым и голубым - типа a+bc, тёмно-синим - вычисление [math]\gamma , v_{1}[/math], жёлтым - умножения (или деления). Обведённая группа операций повторяется независимо n-i раз. Результаты синего, красных и жёлтых кружков, а на последнем шаге и голубого - выходные для алгоритма.

Как уже сказано в описании ядра, основная часть - вычисления на каждом шагу скалярных произведений [math](s,s)[/math] и [math](x,v)[/math] для всех подстолбцов [math]x[/math] справа от текущего, а также векторные операции [math]x'=x-\frac{(x,v)}{\gamma}v[/math]. При этом, однако, строгая последовательность выполнения этих трёх подшагов не обязательна, в силу связи получаемых векторов [math]s[/math] и [math]v[/math] можно одновременно с [math](s,s)[/math] вычислять и произведения [math](x,s)[/math] с последующим выражением через них [math](x,v)[/math]. Это позволяет почти вдвое уменьшать критический путь графа алгоритма.

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

Последовательность выполнения алгоритма обычно записывается как последовательное "обнуление" поддиагональных элементов столбцов, начиная с 1-го столбца и заканчивая предпоследним [math](n-1)[/math]-м.

При этом в каждом "обнуляемом" [math]i[/math]-м столбце "обнуляются" сразу все его поддиагональные элементы одновременно, с [math](i+1)[/math]-го до [math]n[/math]-го.

Каждое "обнуление" [math]i[/math]-го столбца состоит из двух шагов: а) вычисление параметров матрицы отражения [math]U_{i}[/math] такой, чтобы при умножении на неё слева "обнулились" все поддиагональные его элементы; б) умножение слева матрицы отражения [math]U_{i}[/math] на текущую версию матрицы.

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

В последовательной версии основная сложность алгоритма определяется прежде всего вычислениями скалярных произведений векторов, а также модификаций векторов вида [math]x'=x-\alpha v[/math], причем над векторами убывающей по ходу алгоритма размерности. Они, если не учитывать возможную разреженность, составляют (в главном члене) по [math]2n^3/3[/math] операций действительного умножения и сложения/вычитания.

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

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

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

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

Для понимания ресурса параллелизма в разложении матрицы порядка [math]n[/math] методом Хаусхолдера нужно рассмотреть критический путь графа.

Как видно из описания разных вершин, вычисления при "обнулении" [math]i[/math]-го столбца параметров отражения и скалярных произведений состоят из основной части - ветви длиной по [math]n-i[/math] умножений и сложений - и коррекции вычислений, которые составляют [math]O(1)[/math] операций.

Поэтому по грубой (без членов низших порядков) оценке критический путь метода Хаусхолдера будет идти через [math]\frac{n^2}{2}[/math] умножений и [math]\frac{n^2}{2}[/math] сложений/вычитаний.

Поэтому в параллельном варианте, как и в последовательном, основную долю требуемого для выполнения алгоритма времени будут определять операции вида [math]a+bc[/math].

При классификации по высоте ЯПФ, таким образом, метод Хаусхолдера относится к алгоритмам с квадратичной сложностью. При классификации по ширине ЯПФ его сложность будет также квадратичной (без расширения ряда ярусов, связанных с векторными операциями сложения, пришлось бы увеличить вдвое длину критического пути; при таком расширении сложность по ширине ЯПФ станет линейной).

Надо сказать, что здесь в оценках речь идёт именно о классическом способе реализации метода Хаусхолдера. Даже использование схем сдваивания или последовательно-параллельных для вычисления скалярных произведений уменьшает критический путь с квадратичного до степени 3/2 или линейно-логарифмического. Однако все эти широко распространённые методы пока не дали возможности снизить критический путь метода Хаусхолдера до линейного (как, скажем, у метода Гивенса).

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

Входные данные: плотная квадратная матрица [math]A[/math] (элементы [math]a_{ij}[/math]).

Объём входных данных: [math]n^2[/math].

Выходные данные: правая треугольная матрица [math]R[/math] (ненулевые элементы [math]r_{ij}[/math] в последовательном варианте хранятся в элементах исходной матрицы [math]a_{ij}[/math]), унитарная (ортогональная) матрица Q - как произведение матриц Хаусхолдера (отражения) (их вектора нормалей к плоскостям отражения в последовательном варианте хранятся в поддиагональных элементах исходной матрицы [math]a_{ij}[/math] и в одном дополнительном столбце размерности n).

Объём выходных данных: [math]n^2+n[/math].

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

Соотношение последовательной и параллельной сложности, как хорошо видно, является линейным, что даёт определённый стимул для распараллеливания. Однако у наискорейшей ЯПФ ширина квадратична, что указывает на дисбаланс между загруженностями устройств при попытке её реально запрограммировать. Поэтому более практично даже при хорошей (быстрой) вычислительной сети оставить количество устройств (например, узлов кластера) линейным по размеру матрицы, что удвоит критический путь реализуемой ЯПФ.

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

Алгоритм в рамках выбранной версии полностью детерминирован.

Вычислительная погрешность в методе отражений (Хаусхолдера) растет линейно, как и в методе Гивенса (вращений).

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

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

В варианте с кратчайшим критическим путём графа алгоритма (с использованием зависимости между обнуляемым вектором и направляющим вектором отражения) метод Хаусхолдера (отражений) QR-разложения квадратной вещественной матрицы на Фортране 77 можно записать так:

	DO  I = 1, N-1
            DO K = I, N
              SX(K)=A(N,I)*A(N,K)
            END DO
            DO J = N-1, I, -1
               DO K = I, N
                  SX(K)=SX(K)+A(J,I)*A(J,K)
               END DO
            END DO

               ALPHA = SQRT (SX(I))
               IF (A(I,I).NE.0.) THEN
                  BETA = 1./ALPHA
                  DO J = I+1, N
                    A(J,I)=A(J,I)*BETA
                  END DO
                  SX(I) = A(I,I)*BETA+SIGN(1.,A(I,I))                     
                  A(I,I)=ALPHA
                  G=1./ABS(SX(I)) ! 1/gamma
                  DO K = I+1, N
                     SX(K)=SX(K)*BETA*G+SIGN(A(I,K),SX(I))
                     A(I,K) = A(I,K)+SX(K)*SX(I)
                     DO J = I+1, N
                       A (J,K) = A(J,K)-A(J,I)*SX(K)
                     END DO
                  END DO
               ELSE
                  IF (ALPHA.NE.0.) THEN
                     BETA = 1./ALPHA
                     DO J = I+1, N
                       A(J,I)=A(J,I)*BETA
                     END DO
                     SX(I) = -1. 
                     A(I,I)=ALPHA
                     G=1.! 1/gamma
                     DO K = I+1, N
                        SX(K)=SX(K)*BETA*G+SIGN(A(I,K),SX(I))
                        A(I,K) = A(I,K)+SX(K)*SX(I)
                        DO J = I+1, N
                          A (J,K) = A(J,K)-A(J,I)*SX(K)
                        END DO
                     END DO
                  ELSE
                     SX(I)=1
                     G=2.
                     DO K = I+1, N
                        SX(K)=2.
                        A(I,K) = A(I,K)-SX(K)
                     END DO
                  END IF
               END IF               

		
	END DO

В этом варианте R расположена в верхнем правом треугольнике массива A, направляющие вектора матриц отражений размещены в поддиагональных элементах соответствующих столбцов, а их первые элементы - в элементах массива SX.

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

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

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

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

2.2.1.1 Структура обращений в память и качественная оценка локальности
Рисунок 2. Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы, вещественный вариант. Общий профиль обращений в память

На рис. 2 представлен профиль обращений в память для реализации вещественного варианта метода Хаусхолдера (отражений) QR-разложения квадратной матрицы. Из данного рисунка можно сделать несколько выводов. Во-первых, видно, что в среднем к каждому элементу выполняется примерно 90 обращений (отношение общего числа обращений к числу задействованных данных), что является достаточно хорошим показателем в плане повторного использования данных. Далее, явно видна итерационная структура профиля, причем на каждой следующей итерации обращения к небольшой части первых элементов перестают выполняться. Это говорит о том, что ближе к концу профиля задействовано меньше данных, что хорошо сказывается на пространственной локальности. Однако в то же время видно, что на каждой следующей итерации выполняется меньше обращений к одним и тем же локальным областям, то есть временная локальность, по всей вероятности, становится ниже. Рассмотрим подробнее первую итерацию, что более детально оценить их строение.

На рис. 3 представлена одна итерация общего профиля (выделенный зеленый фрагмент на рис. 2). Ее можно условно разбить на 5 фрагментов, каждую из которых рассмотреть отдельно. Фрагменты 2 и 5 устроены просто – элементы перебираются с большим шагом по памяти, при достижении конца данных перебор повторяется, только с небольшим сдвигом. Про такие фрагменты можно сказать, что они обладают низкой пространственной локальностью, поскольку шаг по памяти достаточно большой, и низкой временной локальностью, так как данные повторно не используются. Однако отметим, что данные фрагменты состоят из небольшого числа обращений, что означает, что они не так сильно будут снижать общую локальность профиля.

Рисунок 3. Профиль обращений, одна итерация

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

Рисунок 4. Профиль обращений, одна итерация, фрагмент 3

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

Далее рассмотрим фрагмент 4. Чтобы лучше понять локальную структуру, детально изучим небольшую его часть (рис. 5). Из данного рисунка хорошо видно, что здесь профиль устроен очень просто и представляет собой обычный последовательный перебор с шагом 1 (с небольшими и нечастыми сдвигами, которые не влияют на локальность этого фрагмента). Такой фрагмент обладает высокой пространственной, но низкой временной локальностью.

Рисунок 5. Профиль обращений, одна итерация, часть фрагмента 4

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

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

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

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


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

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

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

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

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

Проведём исследование масштабируемости вширь реализации согласно методике. Исследование проводилось на суперкомпьютере "Ломоносов-2" Суперкомпьютерного комплекса Московского университета.

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

  • число процессоров [16:100] с шагом n^2;
  • размер матрицы [1000 : 16000] с шагом 1000.

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

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

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

Рисунок 7. Параллельная реализация метода Хаусхолдера QR-разложения. Изменение производительности в зависимости от числа процессоров и размера матрицы.
Рисунок 8. Параллельная реализация метода Хаусхолдера QR-разложения. Изменение эффективности в зависимости от числа процессоров и размера матрицы.

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

Для проведения экспериментов использовалась реализация метода Хаусхолдера, представленная в пакете SCALAPACK библиотеки Intel MKL (метод pdgehrd). Все результаты получены на суперкомпьютере «Ломоносов-2». Использовались процессоры Intel Xeon E5-2697v3, сеть FDR Infiniband а также компилятор Intel с опциями, рекомнедованными Intel MKL Advisor. На рисунках показана эффективность реализации метода Хаусхолдера (отражений) QR-разложения квадратной матрицы для размерности матрицы 20000, запуск проводился на 100 ядрах, выполнялось 11 итераций.

Рисунок 9. График загрузки CPU при выполнении метода Хаусхолдера (отражений) QR-разложения квадратной матрицы

На графике загрузки процессора видно, что почти все время работы программы средний уровень загрузки составляет около 85%. Это хорошая картина для программ, запущенных c использованием технологии Hyper Threading т.е. всех виртуальных ядер доступных для данной архитектуры. К концу итерации виден более высокий уровень интенсивности вычислений в срванении со началом и серединой вычислений.

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

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

Рисунок 11. График кэш-промахов L1 в секунду при работе метода Хаусхолдера (отражений) QR-разложения квадратной матрицы

На графике кэш-промахов первого уровня видно, что число промахов достаточно низкое и находится на уровне 5 млн/сек в среднем по всем узлам. Интересен факт снижения числа промахов к концу каждой итерации, хотя снижение и не очень значительное, но вполне заметное.

Рисунок 12. График кэш-промахов L1 в секунду при работе метода Хаусхолдера (отражений) QR-разложения квадратной матрицы

На графике кэш-промахов второго уровня видно, что число промахов достаточно тоже низкое и находится на уровне 1 млн/сек в среднем по всем узлам. На графике промахов второго уровня факт снижения числа промахов к концу итерации проявляется более явно и снижение и более значительное и заметное (с 1.2 млн до 0.5 млн.).

Рисунок 13. График кэш-промахов L3 в секунду при работе метода Хаусхолдера (отражений) QR-разложения квадратной матрицы

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

Рисунок 13. График скорости передачи по сети Infiniband в байт/сек при работе метода Хаусхолдера (отражений) QR-разложения квадратной матрицы

На графике скорости передачи данных по сети Infiniband наблюдается достаточно низкая интенсивность использования коммуникационной сети на каждой итерации. Причем к концу каждой итерации интенсивность передачи данных сильно снижается. Это указывает на бОльшую необходимость в обмене данными между процессами в начале итерации и большую локальность вычислений к концу итерации. Наблюдается так же разбалансированность пересылок между максимальной и минимальной интенсивностью среди использованных узлов.

Рисунок 14. График скорости передачи по сети Infiniband в пакетах/сек при работе метода Хаусхолдера (отражений) QR-разложения квадратной матрицы

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

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

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

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

2.7 Существующие реализации алгоритма

Большинство пакетов от LINPACKа и LAPACKa до SCALAPACKa используют для QR-разложения матриц именно метод Хаусхолдера, правда, в различных модификациях (обычно с использованием BLAS). Существует большая подборка исследовательских работ по блочным версиям.

3 Литература

  1. В.В.Воеводин, Ю.А.Кузнецов. Матрицы и вычисления. М.: Наука, 1984.
  2. Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.