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

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

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

Содержание

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

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

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

На каждом шаге алгоритма две строки преобразуемой матрицы подвергаются преобразованию вращения. При этом параметр такого преобразования подбирается так, чтобы один из поддиагональных элементов преобразуемой матрицы стал нулевым. Сначала в 0 последовательно обращаются элементы 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 Математическое описание алгоритма

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

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

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

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

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

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

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

Рисунок 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 Входные и выходные данные алгоритма

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 Структура обращений в память и качественная оценка локальности
2.2.1.2 Количественная оценка локальности

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

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

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

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

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

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

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

3 Литература

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