Итерация алгоритма dqds: различия между версиями
| [непроверенная версия] | [непроверенная версия] |
| Строка 67: | Строка 67: | ||
Входные и выходные данные: <math>q_{j}, \; j\in [1,n], \; e_{k}, \; k\in [1,n-1] </math> - квадраты элементов главной и побочной диагонали входной матрицы <math>B</math>, <math> \widehat{q}_j , \; \widehat{e}_k </math> - то же для вычисляемой матрицы <math>\widehat{B}.</math>. | Входные и выходные данные: <math>q_{j}, \; j\in [1,n], \; e_{k}, \; k\in [1,n-1] </math> - квадраты элементов главной и побочной диагонали входной матрицы <math>B</math>, <math> \widehat{q}_j , \; \widehat{e}_k </math> - то же для вычисляемой матрицы <math>\widehat{B}.</math>. | ||
| − | Формулы метода выглядят следующим | + | Формулы метода выглядят следующим образом: |
Версия 15:28, 7 июня 2016
| Алгоритм dqds нахождения сингулярных чисел двухдиагональной матрицы | |
| Последовательный алгоритм | |
| Последовательная сложность | [math]5n-4[/math] |
| Объём входных данных | [math]2n[/math] |
| Объём выходных данных | [math]2n[/math] |
| Параллельный алгоритм | |
| Высота ярусно-параллельной формы | [math]4n-3[/math] |
| Ширина ярусно-параллельной формы | [math]2[/math] |
Основные авторы описания: А.Ю.Чернявский
Содержание
- 1 Свойства и структура алгоритмов
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритмов
- 3 Литература
1 Свойства и структура алгоритмов
1.1 Общее описание алгоритма
Итерация алгоритма dqds является одним шагом алгоритма dqds нахождения сингулярных чисел двухдиагональной матрицы.
Сам алгоритм dqds (differential quotient-difference algorithm with shifts)[1][2] строит последовательность двухдиагональных матриц, сходящуюся к диагональной матрице, содержащей квадраты искомых сингулярных чисел. Его особенностью является высокая точность решния задачи. Вычислительным ядром алгоритма является именно внутренняя итерация, вне итераций происходит подбор сдвига [math]\delta[/math] (параметер итерации, см. математическое описание алгоритма), отслеживание сходимости, а также применение различных оптимизационных "хитростей". Отметим, что внеитерационная часть алгоритма не существенна с точки зрения структуры вычислений, т.к. основные затраты ложаться на вычисления внутри итерации. Подробности и варианты внеитерационной части, а также анализ сходимости можно найти в соответствующей литературе [3] [4] [5].
1.2 Математическое описание алгоритма
1.2.1 Вспомогательные свдения
Для понимания математических основ dqds-итерации полезно рассмотреть кратко её вывод, частично отражающий и историю возникновения алгоритма (подробности можно найти в [1]). За основу dqds-алгоритма удобно взять так называемую LR-итерацию, предшествующую хорошо-известной QR-итерации. LR-алгоритм, начиная с входной матрицы [math]T_0\gt 0,[/math] строит сходящуюся последовательность подобных [math]T_0[/math] матриц [math]T_i\gt 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+1}=B_iB_i^T+\tau_i^2I.[/math]
Отметим, что два шага LR-итерации с нулевым сдвигом эквивалентны одному шагу QR-итерации. Итерационная процедура приводит матрицу к диагональному виду, тем самым вычисляя собственные значения исходной матрицы. LR-алгоритм достаточно легко может быть переформулирован с заметными упрощениями для задачи поиска сингулярных значений двухдиагональных матриц. А именно, будем вычислять последовательность двухдиагональных матриц [math]B_i[/math] без непосредственного вычисления [math]T_i[/math](которые в данном случае будут трехдиагональными). Пусть матрица [math]B_i[/math] имеет диагональные элементы [math]a_1 \ldots a_n[/math] и наддиагональные элементы [math]b_1 \ldots b_{n-1}[/math], а матрица [math]B_{i+1}[/math] - диагональные элементы [math]\widehat{a}_1 \ldots \widehat{a}_n[/math] и наддиагональные элементы [math]\widehat{b}_1 \ldots \widehat{b}_{n-1}.[/math] Тогда шаг LR-итерации в терминах матриц [math]B_i[/math] можно привести к простому циклу, пробегающему значения [math]j[/math] от [math]1[/math] до [math]n-1:[/math]
- [math] \widehat{a}^2_j = a^2_j+b^2_j-\widehat{b}^2_{j-1}-\delta [/math]
- [math] \widehat{b}^2_j = b^2_j\dot (a^2_{j+1}/a^2_j) [/math]
и вычислению [math]\widehat{a}^2_n = a^2_n-\widehat{b}^2_{n-1}-\delta.[/math] Очевидно, что работу с извлечением квадратов выгодно вести лишь после окончания работы алгоритма, поэтому можно ввести замену [math]q_j=a^2_j,\; e_j=b^2,[/math] что в итоге приводит к так называемому алгоритму qds. Формулы алгоритма следующие:
- [math] \widehat{q}_j = q_j + e_j - \widehat{e}_{j-1} - \delta, \quad j \in [1,n-1] [/math]
- [math] \widehat{e}_j = e_j \cdot q_{j+1} / \widehat{q}_j, \quad j \in [1,n-1] [/math]
- [math] \widehat{q}_n = q_n - \widehat{e}_{n-1} - \delta. [/math]
Здесь [math]q_j, \; j \in [1,n][/math] и [math]e_j, \; j \in [1,n-1][/math] - квадраты элемнтов главной и верхней побочной диагонали соответственно. Крышка означает выходные переменные, а
[math]\delta[/math] - сдвиг (параметр алгоритма). Такая математическая запись наиболее компактна и соответсвует так называемой qds-итерации.
1.2.2 Математическое описание итерации алгоритма dqds
Представим теперь математическую запись, приближенную к dqds-итерации (с математической точки зрения qds и dqds-итерации эквивалентны) с введенными вспомогательными переменными
[math]t_j[/math] и [math]d_j.[/math] Итерация алгоритма dqds преобразует входную двухдиагональную матрицу [math]B[/math] в выходную [math]\widehat{B}.[/math]
Входные и выходные данные: [math]q_{j}, \; j\in [1,n], \; e_{k}, \; k\in [1,n-1] [/math] - квадраты элементов главной и побочной диагонали входной матрицы [math]B[/math], [math] \widehat{q}_j , \; \widehat{e}_k [/math] - то же для вычисляемой матрицы [math]\widehat{B}.[/math].
Формулы метода выглядят следующим образом:
- [math] d_1 = q_1 - \delta, \; q_n = d_n [/math]
- для[math] \quad j\in [1,n-1]: [/math]
- [math] \widehat{q}_j = d_j + e_j [/math]
- [math] t_j = q_{j+1}/\widehat{q}_j [/math]
- [math] \widehat{e}_j=e_j \cdot t_j [/math]
- [math] d_{j+1} = d \cdot t - \delta [/math]
1.3 Вычислительное ядро алгоритма
Вычислительным ядром алгоритма является производимое n-1 раз вычисление, содержащие по одной операции сложения, вычитания и деления, а также две операции умножения.
1.4 Макроструктура алгоритма
Алгоритм состоит из отдельного вычисления начального значения вспомогательной переменной [math]d[/math] и последующим (n-1)-кратным выполнением повторяющейся последовательности из 5 операций (+,/,*,*,-).
1.5 Схема реализации последовательного алгоритма
Отметим, что выходные данные сразу могут быть записаны на место входных (это учтено в схеме), также для хранения вспомогательных переменных [math]t_j[/math] и [math]d_j[/math] достаточно двух перезаписываемых переменных.
1. Вычисляется начальное значение вспомогательной переменной [math]d = q_1-\delta.[/math]
2. Производится цикл по j от 1 до n-1, состоящий из:
- 2.1 Вычисляется значение [math]q_j = d + e_j;[/math]
- 2.2 Вычисляется значение вспомогательной переменной [math]t = q_{j+1}/q_j;[/math]
- 2.3 Вычисляется значение [math]e_j = e_j \cdot t;[/math]
- 2.4 Вычисляется значение вспомогательной переменной [math]d = d \cdot t - \delta.[/math]
3. Вычисляется [math]q_n = d.[/math]
Легко заметить, что можно представить вычисления в другой форме, например, в виде qds-итерации (см. Математическое описание dqds-итерации), однако, именно dqds реализация вычисления позволяет достичь высокой точности[1].
1.6 Последовательная сложность алгоритма
Для выполнения одной итерации dqds необходимо выполнить:
- [math]n-1[/math] делений,
- [math]2n-2[/math] умножений,
- [math]2n-1[/math] сложений/вычитаний.
Таким образом одна dqds-итерация имеет линейную сложность.
1.7 Информационный граф
1.8 Ресурс параллелизма алгоритма
Как видно из информационного графа алгоритма, на каждом шаге основного цикла возможно лишь параллельное выполнение операции умножения (2.2) и умножения+сложения (2.4). Это позволяет сократить число ярусов на одной итерации цикла c 5 до 4, а общее число ярусов алгоритма с 5n-4 до 4n-3. Ярусы с операциями умножения состоят из двух операций, остальные же из одной.
1.9 Входные и выходные данные алгоритма
Входные данные: Квадраты элементов основной и верхней побочной диагонали двухдиагональной матрицы (вектора [math]q[/math] длины n и [math]e[/math] длины n-1), а также параметр сдвига [math]\delta[/math].
Объём входных данных: [math]2n[/math].
Выходные данные: Квадраты элементов основной и верхней побочной диагонали выходной двухдиагональной матрицы.
Объём выходных данных: [math]2n-1[/math].
1.10 Свойства алгоритма
Соотношение последовательной и параллельной сложности при наличии возможности параллельного выполнения операций умножения составляет [math]\frac{5n-4}{4n-3}[/math], т.е. алгоритм плохо распараллеливается.
Вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных - константа.
Описываемый алгоритм является полностью детерминированным.
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.
- ↑ Fernando K. V., Parlett B. N. Accurate singular values and differential qd algorithms //Numerische Mathematik. – 1994. – Т. 67. – №. 2. – С. 191-229.
- ↑ Parlett B. N., Marques O. A. An implementation of the dqds algorithm (positive case) //Linear Algebra and its Applications. – 2000. – Т. 309. – №. 1. – С. 217-259.
- ↑ Aishima K. et al. On convergence of the DQDS algorithm for singular value computation //SIAM Journal on Matrix Analysis and Applications. – 2008. – Т. 30. – №. 2. – С. 522-537.
