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

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

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


Метод Гивенса (вращений) QR-разложения квадратной матрицы
Последовательный алгоритм
Последовательная сложность [math]2n^3[/math]
Объём входных данных [math]n^2[/math]
Объём выходных данных [math]n^2[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы [math]11n-16[/math]
Ширина ярусно-параллельной формы [math]O(n^2)[/math]

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

Содержание

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

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

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


номера столбцов: [math]\begin{matrix} \ _{i-1}\quad _i\quad _{i+1} & \ & _{j-1}\ \ _j\quad _{j+1}\end{matrix}[/math]

[math]T_{ij} = \begin{bmatrix} 1 & \cdots & 0\quad 0\quad 0 & \cdots & 0\quad 0\quad 0 & \cdots & 0 \\ \vdots & \ddots & \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & 1\quad 0\quad 0 & \cdots & 0\quad 0\quad 0 & \cdots & 0 \\ 0 & \cdots & 0\quad c\quad 0 & \cdots & 0\ -s\ 0 & \cdots & 0 \\ 0 & \cdots & 0\quad 0\quad 1 & \cdots & 0\quad 0\quad 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & \cdots & 0\quad 0\quad 0 & \cdots & 1\quad 0\quad 0 & \cdots & 0 \\ 0 & \cdots & 0\quad s\quad 0 & \cdots & 0\quad c\quad 0 & \cdots & 0 \\ 0 & \cdots & 0\quad 0\quad 0 & \cdots & 0\quad 0\quad 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0\quad 0\quad 0 & \cdots & 0\quad 0\quad 0 & \cdots & 1 \\ \end{bmatrix} [/math]

может быть определена парой индексов и одним параметром. Это позволяет в стандартной реализации метода Гивенса хранить результаты разложения на месте матрицы [math]A[/math] без использования дополнительных массивов. [math]QR[/math]-разложение матрицы [math]A[/math] может быть использовано в разных целях: как для решения СЛАУ вида [math]Ax = b[/math], так и в качестве одного из шагов так называемого [math]QR[/math]-алгоритма нахождения собственных чисел матрицы.

На каждом шаге алгоритма две строки преобразуемой матрицы подвергаются преобразованию вращения. При этом параметр такого преобразования подбирается так, чтобы один из поддиагональных элементов преобразуемой матрицы стал нулевым. Сначала в нуль последовательно обращаются элементы 1-го столбца, затем 2-го, и т.д. до n-1-го, после чего получившаяся матрица - это и есть матрица [math]R[/math]. Сам шаг алгоритма распадается на две части: подбор параметра вращения и выполнение вращения над двумя строками матрицы. Поскольку слева от "рабочего" столбца элементы вращаемых строк равны 0, то там выполнять вращение не нужно. Вращение элементов строк в текущем столбце выполняется одновременно с подбором параметра вращения. Таким образом, вторая часть шага заключается в выполнении вращения двумерных векторов, составленных из элементов вращаемых строк в столбцах справа от "рабочего". Каждое такое вращение эквивалентно перемножению двух комплексных чисел (или в сумме 4 умножениям, 1 сложению и 1 вычитанию действительных), одно из которых имеет модуль 1. Подбор параметра вращения по двум элементам "рабочего" столбца является более сложной операцией, что связано в том числе и с необходимостью минимизации влияния ошибок округления. Обычно для хранения матрицы вращения используется тангенс половинного угла [math]t[/math], с которым простыми формулами (т.н. "боевыми формулами тригонометрии") связаны косинус [math]c[/math] и синус [math]s[/math] самого угла:

[math]c = (1 - t^2)/(1 + t^2), s = 2t/(1 + t^2)[/math]

Именно [math]t[/math] обычно и хранится в соответствующем элементе массива.

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

Для получения [math]QR[/math]-разложения квадратной матрицы [math]A[/math] последнюю приводят к правой треугольной ([math]R[/math] - от слова right) последовательностью умножений её слева на матрицы вращения [math]T_{1 2}, T_{1 3}, ..., T_{1 n}, T_{2 3}, T_{2 4}, ..., T_{2 n}, ... , T_{n-2 n}, T_{n-1 n}[/math].

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

В конце процесса имеем [math]R=T_{n-1 n}T_{n-2 n}T_{n-2 n-1}...T_{1 3}T_{1 2}A[/math].

Так как матрицы вращения унитарны, то естественно получается [math]Q=(T_{n-1 n}T_{n-2 n}T_{n-2 n-1}...T_{1 3}T_{1 2})^* =T_{1 2}^* T_{1 3}^* ...T_{1 n}^* T_{2 3}^* T_{2 4}^* ...T_{2 n}^* ...T_{n-2 n}^* T_{n-1 n}^*[/math] и [math]A=QR[/math].

В вещественном случае матрицы вращения ортогональны и [math]Q=(T_{n-1 n}T_{n-2 n}T_{n-2 n-1}...T_{1 3}T_{1 2})^T =T_{1 2}^T T_{1 3}^T ...T_{1 n}^T T_{2 3}^T T_{2 4}^T ...T_{2 n}^T ...T_{n-2 n}^T T_{n-1 n}^T[/math].

Для завершения математического описания остаётся записать вычисление матрицы вращения [math]T_{i j}[/math] и формулы её умножения на текущую промежуточную матрицу.

Пусть в позиции [math](i,i)[/math] преобразуемой матрицы стоит число [math]x[/math], а в позиции [math](i,j)[/math] - число [math]y[/math].

Тогда для минимизации влияния ошибок округления сначала вычисляется равномерная норма вектора [math]z = max (|x|,|y|)[/math].

Если она равна 0, то поворот не требуется: [math]t=s=0, c=1[/math].

Если [math]z=|x|[/math], то вычисляем [math]y_1=y/x[/math] и далее [math]c = \frac {1}{\sqrt{1+y_1^2}}[/math], [math]s=-c y_1[/math], [math]t=\frac {1-\sqrt{1+y_1^2}}{y_1}[/math], а в позиции [math](i,i)[/math] после поворота будет число [math]x \sqrt{1+y_1^2}[/math].

Если же [math]z=|y|[/math], то вычисляем [math]x_1=x/y[/math] и далее [math]t=x_1 - x_1^2 sign(x_1)[/math], [math]s=\frac{sign(x_1)}{\sqrt{1+x_1^2}}[/math], [math]c = s x_1[/math], а в позиции [math](i,i)[/math] после поворота будет число [math]y \sqrt{1+x_1^2} sign(x)[/math].

После получения элементов [math]c[/math] и [math]s[/math] матрицы вращения [math]T_{i j}[/math] запись преобразования для каждого столбца правее [math]i[/math]-го выглядит просто. Если в [math]k[/math]-м столбце в позиции c номером [math]i[/math] стоит x, а в позиции с номером [math]j[/math] стоит y, то их новые значения будут [math]cx - sy[/math] и [math]sx + cy[/math], соответственно, как если бы комплексное число с действительной частью [math]x[/math] и мнимой [math]y[/math] умножили на комплексное число [math](c,s)[/math].

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

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

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

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

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

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

При этом в рамках "обнуления" i-го столбца последовательно "обнуляются" его элементы, начиная с i+1-го и заканчивая n-м.

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

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

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

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

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

Макрограф алгоритма изображён на рисунке 1, графы макровершин - на последующих.

Рисунок 1. Граф алгоритма без отображения входных и выходных данных. n=4. F1 - операция вычисления параметров поворота, F2 - выполнение поворота 2-мерного вектора (эквивалентно операции умножения двух комплексных чисел).

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

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

Как видно из описания подграфов, макровершина вычисления параметров поворота F1 в графе намного "весомее" вершины поворота F2. А именно, в макровершине поворота критический путь - всего-навсего одно умножение (их 4, но все можно выполнить параллельно) и одно сложение/вычитание (их 2, но они тоже параллельны). А вот макровершина вычисления параметров даёт в худшем случае критический путь , в котором одно вычисление квадратного корня, и по две операции деления, умножения, сложения/вычитания.

Поэтому по грубой оценке критический путь в основном будет идти через [math]2n-3[/math] макровершины вычисления параметров F1, а также через [math]n-1[/math] макровершин поворотов. В сумме это даёт критический путь, проходящий по [math]2n-3[/math] операциям вычисления квадратного корня, [math]4n-6[/math] делениям, по [math]5n-7[/math] умножениям и сложениям/вычитаниям. В макрографе, представленном на рисунке, один из критических путей может быть представлен, например, как прохождение по "верхней линии" вершин F1 (это [math]n-1[/math] штук) и затем попеременное выполнение F2 и F1 ([math]n-2[/math] раза), после чего ещё одно выполнение F2.

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

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

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]t_{ij}[/math] в последовательном варианте хранятся в элементах исходной матрицы [math]a_{ij}[/math]).

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

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

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

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

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

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

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

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

В простейшем варианте метод Гивенса (вращений) QR-разложения квадратной вещественной матрицы на Фортране можно записать так:

	DO  I = 1, N-1
		DO J = I+1, N
                    CALL PARAMS (A(I,I), A(J,I), C, S)
                    DO K = I+1, N
                       CALL ROT2D (C, S, A(I,K), A(J,K))
                    END DO
                END DO
	END DO

Подпрограмма вращения ROT2D может быть записана, если операции с комплексными числами реализованы в трансляторе хорошо, так:

	SUBROUTINE ROT2D (C, S, X, Y)
        COMPLEX Z
        REAL ZZ(2)
        EQUIVALENCE Z, ZZ
        Z = CMPLX(C, S)*CMPLX(X,Y)
        X = ZZ(1)
        Y = ZZ(2)
        RETURN
        END

или же так:

	SUBROUTINE ROT2D (C, S, X, Y)
        ZZ = C*X - S*Y
        Y = S*X + C*Y
        X = ZZ
        RETURN
        END

(с точки зрения графа алгоритма они эквивалентны).

Подпрограмма вычисления параметров вращения может быть такова:

	SUBROUTINE PARAMS (X, Y, C, S)
        Z = MAX (ABS(X), ABS(Y))
        IF (Z.EQ.0.) THEN
C OR (Z.LE.OMEGA) WHERE OMEGA - COMPUTER ZERO
            C = 1.
            S = 0.
        ELSE IF (Z.EQ.ABS(X))
            R = Y/X
            RR = R*R
            RR2 = SQRT(1+RR)
            X = X*RR2
            Y = (1-RR2)/R
            C = 1./RR2
            S = -C*R
        ELSE
            R = X/Y
            RR = R*R
            RR2 = SQRT (1+RR)
            X = SIGN(Y, X)*RR
            Y = R - SIGN(RR,R)
            S = SIGN(1./RR2, R) 
            C = S*R  
        END IF
        RETURN
        END

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

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

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

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

На рис. 6 представлен профиль обращений в память для реализации метода Гивенса (вращений) QR-разложения квадратной матрицы, для вещественного варианта. Данный профиль состоит из обращений к одному двумерному массиву, в котором хранится матрица значений. Из данного графика можно увидеть наличие однотипных итераций, из которых профиль и состоит. i-я итерация затрагивает элементы массива, начиная с (i-1)*k, то есть после каждой итерации первые k элементов (значение k нам пока неизвестно) больше не рассматриваются. Каждая итерация состоит из двух параллельно выполняемых частей – последовательного перебора всех элементов, начиная с (i-1)*k, и активного использования первых i*k элементов.

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

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

Рисунок 7. Выделенный фрагмент общего профиля

Рассмотрим фрагмент на рис. 7 более подробно (рис. 8). Здесь уже видны отдельные обращения, и теперь можно с уверенностью сказать, что каждый участок представляет собой последовательный перебор небольшого числа элементов. При этом в верхней части на каждом новом участке перебираются новые элементы, а в нижней – одни и те же данные. Поэтому в целом подобный фрагмент обладает высокой пространственной локальностью (обе части состоят из последовательных переборов), но средней временной локальностью (верхняя часть – очень низкой, нижняя – очень высокой локальностью).

Рисунок 8. Фрагмент общего профиля, дальнейшее приближение

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Поиски по известным пакетам программ показали удивительное явление[2]. Несмотря на то, что в большинстве старых "последовательных" пакетов как правило наличествует и программа QR-разложения методом Гивенса, в новых, рассчитанных на параллельное исполнение, пакетах QR-разложение осуществляют только с помощью метода Хаусхолдера (отражений), который именно по параллельным свойствам проигрывает методу Гивенса. Мы связываем это с тем, что метод Хаусхолдера (отражений) легко записывается через BLAS-процедуры и не требует переупорядочивания обычного хода циклов (разве что при введении блочности, но это не "косая", а координатная нарезка циклов), а вот для параллельной реализации метода Гивенса (вращений) не избежать перезаписи циклов, поскольку максимально параллельный вариант в нём требует применения cкошенного (в узком понимании этого термина) параллелизма. Поэтому временно описание масштабируемости реализаций данного алгоритма невозможно - пока не у чего исследовать эту масштабируемость.

3 Литература

  1. В.В.Воеводин, Ю.А.Кузнецов. Матрицы и вычисления. М.: Наука, 1984.
  2. Фролов А.В., Антонов А.С., Воеводин Вл.В., Теплов А.М. Сопоставление разных методов решения одной задачи по методике проекта Algowiki // Подана на конференцию "Параллельные вычислительные технологии (ПаВТ'2016)".