Классический метод Якоби (вращений) для симметричных матриц с выбором по всей матрице
Авторы описания: А.С.Галкина (входные и выходные данные, математическое описание алгоритма и др.), И.А.Плахов (ресурс параллелизма алгоритма,последовательная сложность алгоритма и др.), А.В.Фролов (части про программную реализацию, а также общая редакция и правка страницы).
Метод Якоби — итерационный алгоритм для вычисления собственных значений и собственных векторов вещественной симметричной матрицы. Карл Густав Якоб Якоби, в честь которого назван этот метод, предложил его в 1846 году. Однако использоваться метод начал только в 1950-х годах с появлением компьютеров.
Свойства и структура алгоритма
Общее описание алгоритма
Суть алгоритма заключается в том, чтобы для заданной симметрической матрице [math]\displaystyle{ A = A^{(0)} }[/math] построить последовательность ортогонально подобных матриц [math]\displaystyle{ A^{(1)},A^{(2)},\dotsc,A^{(m)} }[/math], сходящуюся к диагональной матрице, на диагонали которой стоят собственные значения [math]\displaystyle{ A }[/math]. Для построения этой последовательности применяется специально подобранная матрица вращения [math]\displaystyle{ J_i }[/math], такая что норма наддиагональной части [math]\displaystyle{ \| A^{(i)} \|_{off} =\sqrt{\sum\limits_{1 \le j \lt k \le n} (a_{jk}^{(i)})^2} }[/math] уменьшается при каждом двустороннем вращении матрицы [math]\displaystyle{ A^{(i+1)}={J_i}^TA^{(i)}J_i }[/math].
Математическое описание алгоритма
Исходные данные: симметрическая матрица [math]\displaystyle{ A }[/math] (элементы [math]\displaystyle{ a_{ij} }[/math]).
Вычисляемые данные: диагональная матрица [math]\displaystyle{ \Lambda }[/math] (элементы [math]\displaystyle{ \lambda_{ij} }[/math]).
Норма наддиагональной части [math]\displaystyle{ \| A^{(i)} \|_{off} =\sqrt{\sum\limits_{1 \le j \lt k \le n} (a_{jk}^{(i)})^2} }[/math] уменьшается при каждом двустороннем вращении матрицы [math]\displaystyle{ A^{(i+1)}={J_i}^TA^{(i)}J_i }[/math].
Это достигается выбором максимального по абсолютной величине элемента матрицы [math]\displaystyle{ A^{(i)} }[/math] и его обнулением[1] в матрице [math]\displaystyle{ A^{(i+1)} }[/math]. Если он расположен в j-й строке и k-м столбце, то
[math]\displaystyle{ j }[/math] [math]\displaystyle{ k }[/math]
[math]\displaystyle{ J_i = \begin{matrix} \\ \\ \\ j \\ \\ k \\ \\ \\ \\ \end{matrix} }[/math] [math]\displaystyle{ \begin{bmatrix} & 1 & & & & & & & & \\ & & 1 & & & & & & & \\ & & & \ddots & & & & & & \\ & & & & \cos(\theta) & & -\sin(\theta) & & & \\ & & & & & \ddots & & & & \\ & & & & \sin(\theta) & & \cos(\theta) & & & \\ & & & & & & & \ddots & & \\ & & & & & & & & 1 & \\ & & & & & & & & & 1 \\ \end{bmatrix} }[/math]
Если обозначить [math]\displaystyle{ s = \sin \theta }[/math] и [math]\displaystyle{ c = \cos \theta }[/math], то матрица [math]\displaystyle{ A^{(i+1)} }[/math] состоит из следующих элементов, отличающихся от элементов [math]\displaystyle{ A^{(i)} }[/math]:
- [math]\displaystyle{ \begin{align} a_{jj}^{(i+1)} &= c^2\, a_{jj}^{(i)} - 2\, s c \,a_{jk}^{(i)} + s^2\, a_{kk}^{(i)} \\ a_{kk}^{(i+1)} &= s^2 \,a_{jj}^{(i)} + 2 s c\, a_{jk}^{(i)} + c^2 \, a_{kk}^{(i)} \\ a_{jk}^{(i+1)} &= a_{kj}^{(i+1)} = (c^2 - s^2 ) \, a_{jk}^{(i)} + s c \, (a_{kk}^{(i)} - a_{jj}^{(i)} ) \\ a_{jm}^{(i+1)} &= a_{mj}^{(i+1)} = c \, a_{jm}^{(i)} - s \, a_{km}^{(i)} & m \ne j,k \\ a_{km}^{(i+1)} &= a_{mk}^{(i+1)} = s \, a_{jm}^{(i)} + c \, a_{km}^{(i)} & m \ne j,k \\ a_{ml}^{(i+1)} &= a_{ml}^{(i)} &m,l \ne j,k \end{align} }[/math]
Можно выбрать [math]\displaystyle{ \theta }[/math] так, чтобы [math]\displaystyle{ a_{jk}^{(i+1)} = 0 }[/math] и [math]\displaystyle{ a_{kj}^{(i+1)} = 0 }[/math]. Отсюда следует равенство
- [math]\displaystyle{ \frac{a_{jj}^{(i)} - a_{kk}^{(i)}}{2 a_{jk}^{(i)}} = \frac{c^2 - s^2}{2sc} = \frac{\cos (2\theta)}{\sin (2\theta)} = \operatorname{ctg}(2\theta) \equiv \tau }[/math].
Если [math]\displaystyle{ a_{jj}^{(i)} = a_{kk}^{(i)} }[/math], то выбирается [math]\displaystyle{ \theta = \frac{\pi}{4} }[/math], в противном случае
вводится [math]\displaystyle{ t = \frac{s}{c} = \operatorname{tg}(\theta) }[/math] и тогда [math]\displaystyle{ t^2 - 2t\tau + 1 = 0 }[/math]. Решение квадратного уравнения даёт [math]\displaystyle{ t = \frac{\operatorname{sign}(\tau)}{|\tau| + \sqrt{1+\tau^2}}, c = \frac{1}{\sqrt{1+t^2}}, s = tc }[/math].
Вычисление останавливается, когда выполняются критерии близости к диагональной матрице. Это малость максимального по абсолютной величине внедиагонального элемента матрицы.
Вычислительное ядро алгоритма
Рассматривая отдельную итерацию, можно считать, что вычислительное ядро составляют множественные вычисления элементов матрицы [math]\displaystyle{ a_{jm}^{(i+1)} = a_{mj}^{(i+1)} }[/math] и [math]\displaystyle{ a_{km}^{(i+1)} = a_{mk}^{(i+1)} }[/math], [math]\displaystyle{ m \ne j,k }[/math] в процессе применения матрицы поворота [math]\displaystyle{ J_i }[/math] к матрице [math]\displaystyle{ A }[/math]:
- [math]\displaystyle{ \begin{align} a_{jm}^{(i+1)} &= a_{mj}^{(i+1)} = c \, a_{jm}^{(i)} - s \, a_{km}^{(i)} & m \ne j,k \\ a_{km}^{(i+1)} &= a_{mk}^{(i+1)} = s \, a_{jm}^{(i)} + c \, a_{km}^{(i)} & m \ne j,k, \end{align} }[/math]
каждое из которых повторяется по [math]\displaystyle{ (n-2) }[/math] раза, а также вычисление элементов [math]\displaystyle{ a_{jj}^{(i+1)} }[/math] и [math]\displaystyle{ a_{kk}^{(i+1)} }[/math]:
- [math]\displaystyle{ \begin{align} a_{jj}^{(i+1)} &= c^2\, a_{jj}^{(i)} - 2\, s c \,a_{jk}^{(i)} + s^2\, a_{kk}^{(i)} \\ a_{kk}^{(i+1)} &= s^2 \,a_{jj}^{(i)} + 2 s c\, a_{jk}^{(i)} + c^2 \, a_{kk}^{(i)} \end{align} }[/math]
Макроструктура алгоритма
Основную часть метода составляет процедура применения вращения к матрице [math]\displaystyle{ A }[/math] (здесь и далее верхние индексы, содержащие номер матрицы, опускаются), которая в дальнейшем будет обозначена как Jakobi-Rotation[math]\displaystyle{ (A,j,k) }[/math].
Эту процедуру, в свою очередь, можно разделить на две логические части:
- Определение угла поворота [math]\displaystyle{ \theta }[/math] по элементам матрицы [math]\displaystyle{ a_{jj} }[/math], [math]\displaystyle{ a_{kk} }[/math] и [math]\displaystyle{ a_{jk} }[/math];
- Поворот матрицы [math]\displaystyle{ A }[/math] (изменяются лишь строки и столбцы, соответствующие индексам [math]\displaystyle{ j }[/math] и [math]\displaystyle{ k }[/math]).
Схема реализации последовательного алгоритма
Алгоритм можно описать следующим образом:
1. Выбрать пару индексов j,k так, чтобы элемент в этой позиции был максимальным по абсолютной величине из всех внедиагональных 2. Обратиться к процедуре Jakobi-Rotation (A,j,k) Если A недостаточно близка к диагональной матрице, перейти к шагу 1.
Процедура Jakobi-Rotation (A,j,k) — это следующий алгоритм:
Если |a(j,k)| достаточно мал, вычисление заканчивается. В противном случае выполняются следующие действия:
1. Если a(j,j)==a(k,k), то угол theta = pi/4 .
В остальных случаях находим параметры tau, t, c, s:
tau = (a(j,j)-a(k,k)/2/a(j,k)
t = sign(tau)/(|tau|+sqrt(1+tau^2))
c = 1/sqrt(1+t^2)
s = ct
2. Выполняется поворот матрицы (изменяются лишь строки и столбцы, соответствующие индексам k и j):
A = R*(j,k,theta)AR(j,k,theta), c = cos (theta), s = sin (theta)
Последовательная сложность алгоритма
Для осуществления одной итерации метода Якоби для матрицы размера [math]\displaystyle{ (n\times n) }[/math] требуется выполнить:
Для вычислительного ядра —
- [math]\displaystyle{ 4n+4 }[/math] умножений,
- [math]\displaystyle{ 2n+2 }[/math] сложений.
Для остальной части алгоритма —
- [math]\displaystyle{ 3 }[/math] умножения,
- [math]\displaystyle{ 3 }[/math] деления,
- [math]\displaystyle{ 5 }[/math] сложений (вычитаний),
- [math]\displaystyle{ 2 }[/math] извлечения квадратного корня.
Умножения и сложения (вычитания) составляют основную часть алгоритма.
Информационный граф
Ресурс параллелизма алгоритма
Для выполнения одной итерации требуется последовательно выполнить следующие действия:
- вычислить [math]\displaystyle{ \tau }[/math] (2 операции сложения, 1 операция деления)
- вычислить [math]\displaystyle{ t }[/math] (1 операция сравнения, 1 операция взятия модуля, 2 операции сложения, 1 операция умножения, 1 операция деления, 1 операция извлечения квадратного корня)
- вычислить [math]\displaystyle{ c }[/math] (1 операция сложения, 1 операция умножения, 1 операция деления, 1 операция извлечения квадратного корня) и [math]\displaystyle{ s }[/math] (1 операция умножения)
После этого выполняется ярус, отвечающий за применение поворота к матрице [math]\displaystyle{ A }[/math] с параллельным выполнением [math]\displaystyle{ 2(n-2) }[/math] операций вычисления элементов матрицы [math]\displaystyle{ a_{jm}^{(i+1)} = a_{mj}^{(i+1)} }[/math] и [math]\displaystyle{ a_{km}^{(i+1)} = a_{mk}^{(i+1)} }[/math], [math]\displaystyle{ m \ne j,k }[/math] , а также 2 операций вычисления элементов [math]\displaystyle{ a_{jj}^{(i+1)} }[/math] и [math]\displaystyle{ a_{kk}^{(i+1)} }[/math]. Наибольшее количество операций содержится в вычислении последних двух элементов (по 6 умножений и 3 сложения).
Входные и выходные данные алгоритма
Входные данные: матрица [math]\displaystyle{ A }[/math] (элементы [math]\displaystyle{ a_{ij} }[/math]). Дополнительные ограничения:
- [math]\displaystyle{ A }[/math] – симметрическая матрица, т. е. [math]\displaystyle{ a_{ij}= a_{ji}, i, j = 1, \ldots, n }[/math].
Объём входных данных: [math]\displaystyle{ \frac{n (n + 1)}{2} }[/math] (в силу симметричности достаточно хранить только диагональ и над/поддиагональные элементы). В разных реализациях эта экономия хранения может быть выполнена разным образом.
Выходные данные: вектор собственных значений [math]\displaystyle{ \Lambda }[/math] (элементы [math]\displaystyle{ \lambda_{ii} }[/math]).
Объём выходных данных: [math]\displaystyle{ n }[/math].
Свойства алгоритма
Метод Якоби является самым медленным из имеющихся алгоритмов вычисления собственных значений симметричной матрицы. Тем не менее, он способен вычислять малые собственные числа и отвечающие им собственные векторы с гораздо большей точностью, чем другие методы [1]. Кроме того, он не предполагает первоначального приведения матрицы к трехдиагональной форме.
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов для метода Якоби является линейным.
Вычислительная мощность алгоритма линейна.
Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности: число итераций зависит от входных данных и порогового значения.
Программная реализация алгоритма
Особенности реализации последовательного алгоритма
При выборе максимального по абсолютной величине внедиагонального элемента матрицы перед первым шагом требуется сравнить [math]\displaystyle{ \frac{n (n - 1)}{2} }[/math] элементов. Если не использовать структурных данных, то каждый раз нужно будет повторять этот выбор, что сделало бы именно его узким местом в программе. Однако при выполнении первого выбора обнуляемого элемента можно завести массив номеров, в которых расположены самые большие по модулю элементы всех "строко-столбцов" (матрица симметричная) матрицы. Поскольку при выполнении преобразования двустороннего вращения полностью изменяются только два "строко-столбца", а в остальных изменяется только по 2 элемента, то можно сразу при вычислении изменяющихся элементов проводить выбор наибольшего из них по абсолютной величине в каждом "строко-столбце", что уберёт из алгоритма это узкое место.
Возможные способы и особенности параллельной реализации алгоритма
При выборе максимального по абсолютной величине внедиагонального элемента матрицы перед первым шагом требуется сравнить [math]\displaystyle{ \frac{n (n - 1)}{2} }[/math] элементов, что составляет при использовании последовательно-параллельной схемы O(n) шагов, а при использовании сдваиваний O(log n) шагов. Если не использовать структурных данных, то каждый раз нужно будет повторять этот выбор, что сделало бы именно его узким местом в программе (так как на само двустороннее вращение нужно конечное количество ярусов). Однако при выполнении первого выбора обнуляемого элемента можно завести массив номеров, в которых расположены самые большие по модулю элементы всех "строко-столбцов" (матрица симметричная) матрицы. Поскольку при выполнении преобразования двустороннего вращения полностью изменяются только два "строко-столбца", а в остальных изменяется только по 2 элемента, то можно сразу при вычислении изменяющихся элементов проводить выбор наибольшего из них по абсолютной величине в каждом "строко-столбце". Это не уберёт из алгоритма узкое место, но всё же локализует его на двух преобразуемых "строко-столбцах". Тем не менее, торможение из-за него останется достаточным, чтобы сделать распараллеливание классического метода Якоби слабым его местом даже в сравнении с другим вариантом. Поэтому классический метод Якоби практически не применяется в библиотеках для параллельных вычислительных систем.
По изложенным выше причинам классический метод Якоби практически не применяется в библиотеках для параллельных вычислительных систем.