Метод Хаусхолдера (отражений) приведения матрицы к двухдиагональной форме
| Приведение матрицы к двудиагональному виду методом Хаусхолдера (отражений) | |
| Последовательный алгоритм | |
| Последовательная сложность | [math]\frac{8 n^3}{3}+O(n^2)[/math] |
| Объём входных данных | [math]n^2[/math] |
| Объём выходных данных | [math]n(n + 2)[/math] |
| Параллельный алгоритм | |
| Высота ярусно-параллельной формы | [math]2n^2+O(n)[/math] |
| Ширина ярусно-параллельной формы | [math]n^2[/math] |
Основные авторы описания: А.В.Фролов
Содержание
- 1 Свойства и структура алгоритма
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 2.1 Особенности реализации последовательного алгоритма
- 2.2 Локальность данных и вычислений
- 2.3 Возможные способы и особенности параллельной реализации алгоритма
- 2.4 Масштабируемость алгоритма и его реализации
- 2.5 Динамические характеристики и эффективность реализации алгоритма
- 2.6 Выводы для классов архитектур
- 2.7 Существующие реализации алгоритма
- 3 Литература
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
Метод Хаусхолдера (в советской математической литературе чаще называется методом отражений) используется для приведения симметричных вещественных матриц к двухдиагональному виду, или, что то же самое, для разложения [math]A=QDU^T[/math] ([math]Q, U[/math] - ортогональные, [math]D[/math] — правая двухдиагональная матрица)[1]. При этом матрицы [math]Q, U[/math] хранятся и используются не в своём явном виде, а в виде произведения матриц отражения[2]. Каждая из матриц отражения может быть определена одним вектором. Это позволяет в классическом исполнении метода отражений хранить результаты разложения на месте матрицы A с использованием двух одномерных дополнительных массивов.
В данной статье рассматривается именно классическое исполнение, в котором не используются приёмы типа сдваивания при вычислениях скалярных произведений.
1.2 Математическое описание алгоритма
В методе Хаусхолдера для выполнения разложения матрицы в произведение двухдиагональной и двух ортогональных используются попеременные умножения слева и справа её текущих модификаций на матрицы Хаусхолдера (отражений).
Матрица отражений (Хаусхолдера) - матрица вида [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]i=n-1[/math]) с помощью преобразования отражения "убираются" ненулевые наддиагональные элементы в [math]i[/math]-м столбце, кроме самого левого из них. Таким образом, после [math]n-1[/math] шагов преобразований получается правая двудиагональная [math]D[/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].
Аналогично для правых матриц отражений [math]U=E-\frac{1}{\gamma}vv^*[/math] [math]v[/math] находится через координаты текущей [math]i[/math]-й строки так:
[math]s[/math] - вектор размерности [math]n-i[/math], составленный из элементов [math]i[/math]-й строки, начиная с [math]i+1[/math]-го элемента.
Если [math](s,s)=0[/math], то [math]v=e_{i+1}[/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+1[/math], [math]v_{j}=u_{j-i+2}[/math] при [math]j\gt i+1[/math], а [math]v_{i+1}=1[/math], если [math]u_{1}=0[/math] и [math]v_{i+1}=\frac{u_{1}}{|u_{1}|}(1+|u_{1}|)[/math] для остальных значений. При этом [math]\gamma = 1+|u_{1}|=|v_{i+1}|[/math].
1.3 Вычислительное ядро алгоритма
Основную часть алгоритма составляют вычисления на каждом шагу скалярных произведений [math](s,s)[/math] и [math](x,v)[/math] для всех подстолбцов [math]x[/math] справа от текущего, а также векторные операции [math]x'=x-\frac{(x,v)}{\gamma}v[/math], и затем вычисления на каждом шагу (кроме шага с номером [math]i=n-1[/math]) скалярных произведений [math](s,s)[/math] и [math](x,v)[/math] для всех подстрок [math]x[/math] снизу от текущей, а также векторные операции [math]x'=x-\frac{(x,v)}{\gamma}v[/math]. Это используется при программировании метода во многих библиотеках для его конструирования из стандартных подпрограмм (например, из BLAS).
1.4 Макроструктура алгоритма
Всё приведение состоит из [math]2n-3[/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]-м, попеременно с последовательным "обнулением" наддиагональных (кроме первой кодиагонали) элементов строк, начиная с 1й строки и заканчивая [math](n-2)[/math]-й
При этом в каждом "обнуляемом" [math]i[/math]-м столбце "обнуляются" сразу все его поддиагональные элементы одновременно, с [math](i+1)[/math]-го до [math]n[/math]-го.
Каждое "обнуление" [math]i[/math]-го столбца состоит из двух шагов: а) вычисление параметров матрицы отражения [math]U_{2i-1}[/math] такой, чтобы при умножении на неё слева "обнулились" все поддиагональные его элементы; б) умножение слева матрицы отражения [math]U_{2i-1}[/math] на текущую версию матрицы.
В каждой "обнуляемой" [math]i[/math]-й строке "обнуляются" сразу все его наддиагональные элементы одновременно, с [math](i+2)[/math]-го до [math]n[/math]-го.
Каждое "обнуление" [math]i[/math]-й строки состоит из двух шагов: а) вычисление параметров матрицы отражения [math]U_{2i}[/math] такой, чтобы при умножении на неё справа "обнулились" все (кроме первого) наддиагональные её элементы; б) умножение справа матрицы отражения [math]U_{2i}[/math] на текущую версию матрицы.
1.6 Последовательная сложность алгоритма
В последовательной версии основная сложность алгоритма определяется прежде всего вычислениями скалярных произведений векторов, а также модификаций векторов вида [math]x'=x-\alpha v[/math], причем над векторами убывающей по ходу алгоритма размерности. Они, если не учитывать возможную разреженность, составляют (в главном члене) по [math]4n^3/3[/math] операций действительного умножения и сложения/вычитания.
При классификации по последовательной сложности, таким образом, метод Хаусхолдера относится к алгоритмам с кубической сложностью.
1.7 Информационный граф
На рисунке 1 изображён граф нечётного ([math](2i-1)[/math]го) полушага в привязке к элементам обрабатываемой матрицы. Граф чётного ([math]2i[/math]го) полушага почти аналогичен, только длина ветвей (нарисована слева) будет меньше на 1 ([math]n-i[/math]), и для привязки к элементам матрицы его надо отразить относительно диагонали матрицы. После этого отражения графы полушагов можно положить друг на друга слоями, при этом правые нижние углы должны быть друг над другом.
1.8 Ресурс параллелизма алгоритма
Для понимания ресурса параллелизма в разложении матрицы порядка [math]n[/math] методом Хаусхолдера нужно рассмотреть критический путь графа.
Как видно из описания разных вершин, вычисления при "обнулении" [math]i[/math]-го столбца параметров отражения и скалярных произведений состоят из основной части - ветви длиной по [math]n-i[/math] умножений и сложений - и коррекции вычислений, которые составляют [math]O(1)[/math] операций. Вычисления при "обнулении" [math]i[/math]-й строки параметров отражения и скалярных произведений состоят из основной части - ветви длиной по [math]n-i-1[/math] умножений и сложений - и коррекции вычислений, которые составляют [math]O(1)[/math] операций.
Поэтому по грубой (без членов низших порядков) оценке критический путь метода Хаусхолдера будет идти через [math]n^2[/math] умножений и [math]n^2[/math] сложений/вычитаний.
Поэтому в параллельном варианте, как и в последовательном, основную долю требуемого для выполнения алгоритма времени будут определять операции вида [math]a+bc[/math].
При классификации по высоте ЯПФ, таким образом, метод Хаусхолдера относится к алгоритмам с квадратичной сложностью. При классификации по ширине ЯПФ его сложность будет также квадратичной (без расширения ряда ярусов, связанных с векторными операциями сложения, пришлось бы увеличить вдвое длину критического пути; при таком расширении сложность по ширине ЯПФ станет линейной).
Надо сказать, что здесь в оценках речь идёт именно о классическом способе реализации метода Хаусхолдера. Даже использование схем сдваивания или последовательно-параллельных для вычисления скалярных произведений уменьшает критический путь с квадратичного до степени 3/2 или линейно-логарифмического.
1.9 Входные и выходные данные алгоритма
Входные данные: плотная квадратная матрица [math]A[/math] (элементы [math]a_{ij}[/math]).
Объём входных данных: [math]n^2[/math].
Объём выходных данных: [math]n^2+2n[/math].
1.10 Свойства алгоритма
Соотношение последовательной и параллельной сложности, как хорошо видно, является линейным, что даёт определённый стимул для распараллеливания. Однако у наискорейшей ЯПФ ширина квадратична, что указывает на дисбаланс между загруженностями устройств при попытке её реально запрограммировать. Поэтому более практично даже при хорошей (быстрой) вычислительной сети оставить количество устройств (например, узлов кластера) линейным по размеру матрицы, что удвоит критический путь реализуемой ЯПФ.
При этом вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных, линейна.
Алгоритм в рамках выбранной версии полностью детерминирован.
Вычислительная погрешность в методе отражений (Хаусхолдера) растет линейно, как и в методе Гивенса (вращений).
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
В варианте с кратчайшим критическим путём графа алгоритма (с использованием зависимости между обнуляемым вектором и направляющим вектором отражения) метод Хаусхолдера (отражений) приведения квадратной вещественной матрицы к двудиагональной форме на Фортране 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
IF (I.LT.N-1) THEN
DO K = I, N
SL(K)=A(I,N)*A(K,N)
END DO
DO J = N-1, I+1, -1
DO K = I, N
SL(K)=SL(K)+A(I,J)*A(K,J)
END DO
END DO
ALPHA = SQRT (SL(I))
IF (A(I,I+1).NE.0.) THEN
BETA = 1./ALPHA
DO J = I+2, N
A(I,J)=A(I,J)*BETA
END DO
SL(I) = A(I,I+1)*BETA+SIGN(1.,A(I,I+1))
A(I,I+1)=ALPHA
G=1./ABS(SL(I)) ! 1/gamma
DO K = I+1, N
SL(K)=SL(K)*BETA*G+SIGN(A(K,I+1),SL(I))
A(K,I+1) = A(K,I+1)+SL(K)*SX(I)
DO J = I+2, N
A (K,J) = A(K,J)-A(I,J)*SL(K)
END DO
END DO
ELSE
IF (ALPHA.NE.0.) THEN
BETA = 1./ALPHA
DO J = I+2, N
A(I,J)=A(I,J)*BETA
END DO
SL(I) = -1.
A(I,I+1)=ALPHA
G=1. ! 1/gamma
DO K = I+1, N
SL(K)=SL(K)*BETA*G+SIGN(A(K,I+1),SL(I))
A(K,I+1) = A(K,I+1)+SL(K)*SL(I)
DO J = I+2, N
A (K,J) = A(K,J)-A(I,J)*SL(K)
END DO
END DO
ELSE
SL(I)=1
G=2.
DO K = I+1, N
SL(K)=2.
A(K,I+1) = A(K,I+1)-SL(K)
END DO
END IF
END IF
END IF
END DO
В этом варианте D расположена в диагонали и верхней кодиагонали массива A, направляющие вектора матриц отражений размещены в поддиагональных элементах соответствующих столбцов или надддиагональных - строк, а их первые элементы - в элементах массивов SX и SL.
Обычно же в последовательных версиях коэффициенты модификаций столбцов вычисляются целиком через скалярные произведения после вычислений параметров матрицы отражения. При этом схема чуть проще. Удлиняется критический путь графа, но для последовательных реализаций это неважно.
2.2 Локальность данных и вычислений
К сожалению, в графе, как видно по рисунку, в наличии пучки рассылок, в них неизбежно часть дуг остаются длинными, что отрицательно влияет на локальность вычислений по пространству.
2.2.1 Локальность реализации алгоритма
2.2.1.1 Структура обращений в память и качественная оценка локальности
На рис. 1 представлен профиль обращений в память для реализации метода Хаусхолдера (отражений) приведения матрицы к двухдиагональной форме. Данный профиль обладает явной итерационной структурой, причем видно, что на каждой следующей итерации из рассмотрения отбрасывается небольшая часть элементов из начала массива данных. Также можно заметить, что основная часть данных используется достаточно равномерно в рамках одной итерации, однако в нижней и верхней частях итерации видны небольшие блоки с очень высокой локальностью обращений. Итерации явно устроены по одному принципу, поэтому стоит рассмотреть одну из них более детально.
На рис. 2 показаны обращения в рамках первой итерации общего профиля (выделена на рис. 1 зеленым). Эти обращения можно условно разбить на несколько фрагментов. Фрагменты 1, 2 и 6 практически идентичны (различие заключается только в перестановке местами двух частей этих фрагментов), рассмотрим один из них.
На рис. 3 показан отдельно фрагмент 6, выделенный на рис. 2. Теперь можно определить структуру каждой из двух частей более точно. Первая часть представляет собой последовательный перебор элементов в обратном порядке, причем к каждому элементу выполняется несколько обращений подряд. Во второй части также наблюдается последовательный перебор, однако с другими свойствами: элементы перебираются в прямом порядке, к каждому элементу выполняется только 1 обращение, однако данный перебор повторяется несколько раз подряд, причем каждый раз задействованы одни и те же элементы. Все это позволяет говорить о высокой пространственной локальности фрагмента (поскольку шаг по памяти каждый раз небольшой) и достаточно высокой временной локальности (поскольку данные часто используются повторно).
Свойства локальности во фрагментах 3 и 5 можно оценить и без более детального рассмотрения. Данные в обоих фрагментах перебираются с достаточно большим шагом по памяти, при этом повторно они почти не используются, что говорит о низкой локальности в обоих случаях (хотя надо отметить, что локальность фрагмента 5 немного выше из-за более высокой пространственной локальности).
Теперь перейдем к рассмотрению фрагмента 4, расположенного в центре рис. 2. Поскольку он обладает регулярной структурой, достаточно более детально изучить небольшой фрагмент (представлен на рис. 4, выделен на рис. 2 зеленым). Здесь можно увидеть, что первая часть фрагмента представляет собой практически стандартный последовательный перебор (с небольшими нечастыми сдвигами, которые не влияют на локальность). Такой набор обращений обладает высокой пространственной и низкой временной локальностью. Вторая часть данного фрагмента устроена несколько сложнее – она состоит из L-образных наборов обращений. Одна его область (часть 1 на рис. 3) по строению очень похожа на первую часть; вторая (область 2 на рис. 3) представляет собой набор обращений к одному и тому же элементу, что достаточно сильно повышает локальность данного фрагмента. Из этого можно сделать вывод, что фрагмент 4 обладает высокой пространственной и средней временной локальностью.
Суммируя выше сказанное, можно сказать, что фрагменты 1, 2 и 6 обладают высокой локальностью; фрагменты 3 и 5 – низкой локальностью; фрагмент 4 – высокой пространственной и средней временной локальностью. Можно предположить, что в итоге локальность общего профиля будет не очень высокой, в основном из-за более низкой временной локальности. Однако такое достаточно сложное строение общего профиля затрудняет оценку локальности общего профиля на основе визуального анализа.
2.2.1.2 Количественная оценка локальности
Оценка выполняется на основе характеристики daps, которая оценивает число выполненных обращений (чтений и записей) в память в секунду. Данная характеристика является аналогом оценки flops применительно к работе с памятью и является в большей степени оценкой производительности взаимодействия с памятью, чем оценкой локальности. Однако она служит хорошим источником информации, в том числе для сравнения с результатами по следующей характеристике cvg. На рисунке 5 приведены значения daps для реализаций распространенных алгоритмов, отсортированные по возрастанию (чем больше daps, тем в общем случае выше производительность). Можно увидеть, что оценка производительности в данном случае не очень высока, что подтверждает сделанные ранее выводы. В частности, значение daps для этого примера показывает более низкую производительность по сравнению с реализациями других методов Хаусхолдера (householder_qr householder_qtq) и находится примерно на уровне оценки для реализации метода Холецкого.
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
2.4.1 Масштабируемость алгоритма
2.4.2 Масштабируемость реализации алгоритма
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
В сравнении с методом Гивенса, который имеет естественное двумерное блочное разбиение на основе точечного метода, метод Хаусхолдера из-за худших характеристик локальности (наличие пучков рассылок) и меньшего количества независимых обобщённых развёрток графа не так хорош для реализаций на системах с распределённой памятью, как для систем с общей памятью. Поэтому на массово параллельных системах с распределённой памятью следует применять метод Хаусхолдера (если уж именно его нужно реализовать) не в точечной версии, а в разрабатываемых исследователями блочных вариантах. Следует отметить, что эти варианты - не блочная нарезка описанного метода, а самостоятельные методы. Особенно их применение рекомендуется в случаях с большой разрежённостью матрицы.
2.7 Существующие реализации алгоритма
Большинство пакетов от LINPACKа и LAPACKa до SCALAPACKa используют для приведения матриц к двудиагональному виду именно метод Хаусхолдера, правда, в различных модификациях (обычно с использованием BLAS). Существует большая подборка исследовательских работ по блочным версиям.
