Итерация алгоритма dqds

Алгоритм dqds нахождения сингулярных чисел двухдиагональной матрицы
Последовательный алгоритм
Последовательная сложность	[math]\displaystyle{ 5n-4 }[/math]
Объём входных данных	[math]\displaystyle{ 2n }[/math]
Объём выходных данных	[math]\displaystyle{ 2n }[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]\displaystyle{ 4n-3 }[/math]
Ширина ярусно-параллельной формы	[math]\displaystyle{ 2 }[/math]

Основные авторы описания: А.Ю.Чернявский

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

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

Итерация алгоритма dqds является одним шагом итерационного алгоритма для нахождения сингулярных чисел двухдиагональной матрицы. Основные вычисления алгоритма dqds приходятся именно на внутреннюю итерацию. Внешняя часть алгоритма отвечает за контроль сходимости и выбор сдвигов [math]\displaystyle{ \delta, }[/math] что основано на довольно сложной теории и имеет различные реализации, но в то же время не влияет на вычислительные аспекты. В связи с этим dqds-итерация и dqds-алгоритм рассматриваются отдельно.

Для понимания dqds-итерации полезно рассмотреть кратко её вывод, относительно точно отражающий и историю возникновения алгоритма (подробности можно найти в ). За основу dqds-алгоритма удобно взять так называемую LR-итерацию, предшествующую хорошо-известной QR-итерации. LR-алгоритм, начиная с входной матрицы [math]\displaystyle{ T_0\gt 0, }[/math] строит сходящуюся последовательность подобных [math]\displaystyle{ T_0 }[/math] матриц [math]\displaystyle{ T_i\gt 0, }[/math] итерационно используя следующие три шага:

1. Выбрать сдвиг [math]\displaystyle{ \tau_i }[/math] меньший младшего собственного значения [math]\displaystyle{ T_i. }[/math]
2. Вычислить разложение Холецкого [math]\displaystyle{ T_i-\tau^2_iI=B_i^TB_i, }[/math] где [math]\displaystyle{ B_i }[/math] - верхняя треугольная матрица с положительной диагональю.
3. [math]\displaystyle{ T_{i+1}=B_iB_i^T+\tau_i^2I. }[/math]

Отметим, что два шага LR-итерации эквивалентны одному шагу QR-итерации приводит матрицу к диагональному виду, тем самым вычисляя собственные значения исходной матрицы. Данный алгоритм достаточно легко может быть переформулирован с заметными упрощениями для поиска сингулярных значений двухдиагональных матриц. А именно, будем вычислять последовательность двухдиагональных матриц [math]\displaystyle{ B_i }[/math] без непосредственного вычисления [math]\displaystyle{ T_i. }[/math] Пусть матрица [math]\displaystyle{ B_i }[/math] имеет диагональные элементы [math]\displaystyle{ a_1 \ldots a_n }[/math] и наддиагональные элементы [math]\displaystyle{ b_1 \ldots b_n }[/math], а матрица [math]\displaystyle{ B_{i+1} }[/math] - диагональные элементы [math]\displaystyle{ \widehat{a}_1 \ldots \widehat{a}_n }[/math] и наддиагональные элементы [math]\displaystyle{ \widehat{b}_1 \ldots \widehat{b}_n. }[/math] Тогда шаг LR-итерации в терминах матриц [math]\displaystyle{ B_i }[/math] можно привести к простому циклу, пробегающему значения [math]\displaystyle{ j }[/math] от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ n-1: }[/math]

[math]\displaystyle{ \widehat{a}^2_j = a^2_j+b^2_j-\widehat^2{b}_{j-1}-\delta }[/math]

[math]\displaystyle{ \widehat{b}^2_j = b^2_j\dot (a^2_{j+1}/a^2_j) }[/math]

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

Формулы алгоритма следующие:

[math]\displaystyle{ \widehat{q}_j = q_j + e_j - \widehat{e}_{j-1} - \delta, \quad j \in [1,n-1] }[/math]

[math]\displaystyle{ \widehat{e}_j = e_j \cdot q_{j+1} / \widehat{q}_j, \quad j \in [1,n-1] }[/math]

[math]\displaystyle{ \widehat{q}_n = q_n - \widehat{e}_{n-1} - \delta. }[/math]

Здесь [math]\displaystyle{ q_j, \; j \in [1,n] }[/math] и [math]\displaystyle{ e_j, \; j \in [1,n-1] }[/math] - квадраты элемнтов главной и верхней побочной диагональ соответственно. Крышка означает выходные переменные, а

[math]\displaystyle{ \delta }[/math] - сдвиг (параметр алгоритма). Такая математическая запись наиболее компактна и соответсвует так называемой qds-итерации.

Представим также математическую запись, приближенную к dqds-итерации (с математической точки зрения qds и dqds-итерации эквивалентны) с введенными вспомогательными переменными

[math]\displaystyle{ t_j }[/math] и [math]\displaystyle{ d_j: }[/math]

[math]\displaystyle{ d_1 = q_1 - \delta, \; q_n = d_n }[/math]

для[math]\displaystyle{ \quad j\in [1,n-1]: }[/math]

[math]\displaystyle{ \widehat{q}_j = d_j + e_j }[/math]

[math]\displaystyle{ t_j = q_{j+1}/\widehat{q}_j }[/math]

[math]\displaystyle{ \widehat{e}_j=e_j \cdot t_j }[/math]

[math]\displaystyle{ d_{j+1} = d \cdot t - \delta }[/math]

Отметим, что при [math]\displaystyle{ \delta=0 }[/math] две итерации dqds эквиваленты одной QR-итерации.

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

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

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

Алгоритм состоит из отдельного вычисления начального значения вспомогательной переменной [math]\displaystyle{ d }[/math] и последующим (n-1)-кратным выполнением повторяющейся последовательности из 5 операций (+,/,*,*,-).

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

Отметим, что выходные данные сразу могут быть записаны на место входных (это учтено в схеме), также для хранения вспомогательных переменных [math]\displaystyle{ t_j }[/math] и [math]\displaystyle{ d_j }[/math] достаточно двух перезаписываемых переменных.

1. Вычисляется начальное значение вспомогательной переменной [math]\displaystyle{ d = q_1-\delta. }[/math]

2. Производится цикл по j от 1 до n-1, состоящий из:

2.1 Вычисляется значение [math]\displaystyle{ q_j = d + e_j; }[/math]

2.2 Вычисляется значение вспомогательной переменной [math]\displaystyle{ t = q_{j+1}/q_j; }[/math]

2.3 Вычисляется значение [math]\displaystyle{ e_j = e_j \cdot t; }[/math]

2.4 Вычисляется значение вспомогательной переменной [math]\displaystyle{ d = d \cdot t - \delta. }[/math]

3. Вычисляется [math]\displaystyle{ q_n = d. }[/math]

Легко заметить, что можно представить вычисления в другой форме, например, в виде qds-итерации (см. Математическое описание dqds-итерации), однако, именно dqds реализация вычисления позволяет достичь высокой точности.

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

Для выполнения одной итерации dqds необходимо выполнить:

• [math]\displaystyle{ n-1 }[/math] делений,
• [math]\displaystyle{ 2n-2 }[/math] умножений,
• [math]\displaystyle{ 2n-1 }[/math] сложений/вычитаний.

Таким образом одна dqds-итерация имеет линейную сложность.

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

Рисунок 1. Граф алгоритма для n=4 без отображения входных и выходных данных.

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

Как видно из информационного графа алгоритма, на каждом шаге основного цикла возможно лишь параллельное выполнение операции умножения (2.2) и умножения+сложения (2.4). Это позволяет сократить число ярусов на одной итерации цикла c 5 до 4, а общее число ярусов алгоритма с 5n-4 до 4n-3. Ярусы с операциями умножения состоят из двух операций, остальные же из одной.

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

Входные данные: Квадраты элементов основной и верхней побочной диагонали двухдиагональной матрицы (вектора [math]\displaystyle{ q }[/math] длины n и [math]\displaystyle{ e }[/math] длины n-1), а также параметр сдвига [math]\displaystyle{ \delta }[/math].

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

Выходные данные: Квадраты элементов основной и верхней побочной диагонали выходной двухдиагональной матрицы.

Объём выходных данных: [math]\displaystyle{ 2n-1 }[/math].

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

Соотношение последовательной и параллельной сложности при наличии возможности параллельного выполнения операций умножения составляет [math]\displaystyle{ \frac{5n-4}{4n-3} }[/math], т.е. алгоритм плохо распараллеливается.

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

Описываемый алгоритм является полностью детерминированным.

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

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

Алгоритм на языке Matlab может быть записать так:

d = q(1)-delta;
for j = 1:n-1
    q(j)=d+e(j);
    t=q(j+1)/q(j);
    e(j) = e(j)*t;
    d = d*t-delta;
end
q(n) = d;

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

Легко видеть, что локальность данных высока. Её легко повысить, размещая рядом соответствующие элементы массивов e и q, однако, это не оказывает существенного влияния на производительность в силу константного количества операций относительно объема обрабатываемых данных.

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

Итерация dqds практически полностью последовательна. Единственная возможность - одновременное выполнение операции умножения (2.3) и опрерации (2.4) умножения и сложения, что дает небольшой выигрыш в производительности.

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

Алгоритм не является масштабируемым, максимального эффекта ускорения можно добиться на двух независимых процессорах.

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

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

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

Перечень реализаций приведен на странице алгоритма dqds.