Итерация алгоритма dqds: различия между версиями

[непроверенная версия]

Версия 14:03, 2 июня 2016

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

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

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

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

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

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

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

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

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

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

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

$\widehat{q}_n = q_n - \widehat{e}_{n-1} - \delta.$

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

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

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

$t_j$ и $d_j:$

$d_1 = q_1 - \delta, \; q_n = d_n$

для

$\quad j\in [1,n-1]:$

$\widehat{q}_j = d_j + e_j$

$t_j = q_{j+1}/\widehat{q}_j$

$\widehat{e}_j=e_j \cdot t_j$

$d_{j+1} = d \cdot t - \delta$

Отметим, что при $\delta=0$ две итерации dqds эквиваленты одной QR-итерации.

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

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

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

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

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

Отметим, что выходные данные сразу могут быть записаны на место входных (это учтено в схеме), также для хранения вспомогательных переменных $t_j$ и $d_j$ достаточно двух перезаписываемых переменных.

1. Вычисляется начальное значение вспомогательной переменной $d = q_1-\delta.$

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

2.1 Вычисляется значение

$q_j = d + e_j;$

2.2 Вычисляется значение вспомогательной переменной

$t = q_{j+1}/q_j;$

2.3 Вычисляется значение

$e_j = e_j \cdot t;$

2.4 Вычисляется значение вспомогательной переменной

$d = d \cdot t - \delta.$

3. Вычисляется $q_n = d.$

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

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

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

$n-1$ делений,
$2n-2$ умножений,
$2n-1$ сложений/вычитаний.

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

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

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

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

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

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

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

Объём входных данных: $2n$ .

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

Объём выходных данных: $2n-1$ .

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

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

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

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

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

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

Алгоритм на языке 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;

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

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

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

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

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

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

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

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

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

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

3 Литература

↑ ^{Перейти обратно: 1,0} ^1,1 ^1,2 Деммель Д. Вычислительная линейная алгебра. – М : Мир, 2001.
↑ Hogben L. (ed.). Handbook of linear algebra. – CRC Press, 2006.

[vla-1] {Перейти обратно: 1,0} ^1,1 ^1,2 Деммель Д. Вычислительная линейная алгебра. – М : Мир, 2001.

[hola-2] Hogben L. (ed.). Handbook of linear algebra. – CRC Press, 2006.

[1]

[2]

@@ Строка 21: / Строка 21: @@
 За основу dqds-алгоритма удобно взять так называемую LR-итерацию, предшествующую хорошо-известной QR-итерации. LR-алгоритм, начиная с входной матрицы <math>T_0>0,</math> строит сходящуюся последовательность подобных <math>T_0</math> матриц <math>T_i>0,</math> итерационно используя следующие три шага:
 #Выбрать сдвиг <math>\tau_i</math> меньший младшего собственного значения <math>T_i.</math>
-#Вычислить разложение Холецкого <math>T_i-\tau^2_iI=B_i^TB_i,</math> где <math>B_i</math> -- верхняя треугольная матрица с положительной диагональю.
+#Вычислить разложение Холецкого <math>T_i-\tau^2_iI=B_i^TB_i,</math> где <math>B_i</math> - верхняя треугольная матрица с положительной диагональю.
 #<math>T_{i+1}=B_iB_i^T+\tau_i^2I.</math>