One step of the dqds algorithm
| Алгоритм dqds нахождения сингулярных чисел двухдиагональной матрицы | |
| Sequential algorithm | |
| Serial complexity | [math]5n-4[/math] |
| Input data | [math]2n[/math] |
| Output data | [math]2n[/math] |
| Parallel algorithm | |
| Parallel form height | [math]4n-3[/math] |
| Parallel form width | [math]2[/math] |
Primary authors of this description: A.Yu.Chernyavskiy
Contents
- 1 Properties and structure of the algorithm
- 1.1 General description of the algorithm
- 1.2 Mathematical description of the algorithm
- 1.3 Computational kernel of the algorithm
- 1.4 Macro structure of the algorithm
- 1.5 Implementation scheme of the serial algorithm
- 1.6 Serial complexity of the algorithm
- 1.7 Information graph
- 1.8 Parallelization resource of the algorithm
- 1.9 Input and output data of the algorithm
- 1.10 Properties of the algorithm
- 2 Software implementation of the algorithm
- 2.1 Implementation peculiarities of the serial algorithm
- 2.2 Locality of data and computations
- 2.3 Possible methods and considerations for parallel implementation of the algorithm
- 2.4 Scalability of the algorithm and its implementations
- 2.5 Dynamic characteristics and efficiency of the algorithm implementation
- 2.6 Conclusions for different classes of computer architecture
- 2.7 Existing implementations of the algorithm
- 3 References
1 Properties and structure of the algorithm
1.1 General description of the algorithm
The dqds algorithm iteration implements one step of the dqds algorithm for calculating the singular values of a bi-diagonal matrix.
The full dqds algorithm (differential quotient-difference algorithm with shifts)[1][2] constructs a sequence of bi-diagonal matrices that converges to a diagonal matrix containing the squares of the required singular values. It has the merit of calculating the singular values to high accuracy.
The computational kernel of the algorithm is its inner iteration. Outside of the iterations, shifts [math]\delta[/math] (iteration parameters, see mathematical description of the algorithm) are chosen, convergence is looked after, and various optimization tricks are applied. Note that the non-iterative part of the algorithm is not important with respect to the structure of calculations because the main computational costs fall on the calculations within the iterations. Details concerning the non-iterative part, its variants, and the analysis of convergence can be found in the related literature [3] [4] [5].
1.2 Mathematical description of the algorithm
1.2.1 Auxiliary information
In order to understand the mathematical background of the dqds-iteration, it is helpful to briefly discuss its derivation, which, to same extent, reflects the history of the algorithm (the details can be found in [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 Computational kernel of the algorithm
Вычислительным ядром алгоритма является последовательный расчёт квадратов диагональных ([math]\widehat{q}_j[/math]) и внедиагональных ([math]\widehat{e}_k[/math]) элементов выходной матрицы. Учитывая использование вспомогательных переменных расчёт каждой новой пары содержит по одной операции сложения, вычитания и деления, а также две операции умножения.
1.4 Macro structure of the algorithm
Алгоритм состоит из отдельного вычисления начального значения вспомогательной переменной [math]d,[/math] последующим (n-1)-кратным выполнением повторяющейся последовательности из 5 операций (+,/,*,*,-) для вычисления квадратов диагональных ([math]\widehat{q}_j[/math]) и внедиагональных ([math]\widehat{e}_k[/math]) элементов выходной матрицы и завершающего вычислением крайнего значения [math]\widehat{q}_n[/math].
1.5 Implementation scheme of the serial algorithm
Отметим, что выходные данные сразу могут быть записаны на место входных (это учтено в схеме), также для хранения вспомогательных переменных [math]t_j[/math] и [math]d_j[/math] достаточно двух перезаписываемых переменных. Таким образом элементы главной ([math]q_j[/math]) и побочной ([math]e_k[/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 Serial complexity of the algorithm
Для выполнения одной итерации dqds необходимо выполнить:
- [math]n-1[/math] делений,
- [math]2n-2[/math] умножений,
- [math]2n-1[/math] сложений/вычитаний.
Таким образом одна dqds-итерация имеет линейную сложность.
1.7 Information graph
1.8 Parallelization resource of the algorithm
Как видно из информационного графа алгоритма, на каждом шаге основного цикла возможно лишь параллельное выполнение операции умножения (2.2) и умножения+сложения (2.4). Это позволяет сократить число ярусов на одной итерации цикла c 5 до 4, а общее число ярусов алгоритма с 5n-4 до 4n-3. Ярусы с операциями умножения состоят из двух операций, остальные же из одной.
1.9 Input and output data of the algorithm
Входные данные: Квадраты элементов основной и верхней побочной диагонали двухдиагональной матрицы (вектора [math]q[/math] длины n и [math]e[/math] длины n-1), а также параметр сдвига [math]\delta[/math].
Объём входных данных: [math]2n[/math].
Выходные данные: Квадраты элементов основной и верхней побочной диагонали выходной двухдиагональной матрицы.
Объём выходных данных: [math]2n-1[/math].
1.10 Properties of the algorithm
Соотношение последовательной и параллельной сложности при наличии возможности параллельного выполнения операций умножения составляет [math]\frac{5n-4}{4n-3}[/math], т.е. алгоритм плохо распараллеливается.
Вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных - константа.
Описываемый алгоритм является полностью детерминированным.
2 Software implementation of the algorithm
2.1 Implementation peculiarities of the serial algorithm
Алгоритм на языке 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;
Как говорилось в cхеме реализации последовательного алгоритма, вычисляемые данные записываются сразу на место входных.
2.2 Locality of data and computations
Легко видеть, что локальность данных высока. Её легко повысить, размещая рядом соответствующие элементы массивов e и q, однако, это не оказывает существенного влияния на производительность в силу константного количества операций относительно объема обрабатываемых данных.
2.2.1 Locality of implementation
2.2.1.1 Structure of memory access and a qualitative estimation of locality
На рис. 2 представлен профиль обращений в память для реализации итерации алгоритма dqds. Данный профиль внешне устроен очень просто и состоит из двух параллельно выполняемых последовательных переборов. Отметим, что общее число обращений в память всего в 3 раза больше числа задействованных данных, что говорит о том, что данные повторно используют редко. Зачастую это сигнализирует о достаточно невысокой локальности.
Рассмотрим более детально строение одного из переборов (второй устроен практически идентично). На рис. 3 показан выделенный на рис. 2 небольшой фрагмент 1. Видно, что на каждом шаге в данном переборе выполняется три обращения в память; подобное строение говорит о высокой пространственной локальности, однако низкой временной, поскольку данные повторно используются только по 3 раза.
Поскольку данная структура обращений и составляет общий профиль, то же самое можно говорить и о локальности всего профиля.
2.2.1.2 Quantitative estimation of locality
Оценка выполняется на основе характеристики daps, которая оценивает число выполненных обращений (чтений и записей) в память в секунду. Данная характеристика является аналогом оценки flops применительно к работе с памятью и является в большей степени оценкой производительности взаимодействия с памятью, чем оценкой локальности. Однако она служит хорошим источником информации, в том числе для сравнения с результатами по следующей характеристике cvg.
На рисунке 4 приведены значения daps для реализаций распространенных алгоритмов, отсортированные по возрастанию (чем больше daps, тем в общем случае выше производительность). Результат получен достаточно неожиданный – производительность работы с памятью очень невелика. Оценка daps для данного профиля сравнима с реализациями самых неэффективных алгоритмов в части работы с памятью – тестов RandomAccess и неэффективных вариантов обычного перемножения матриц. Отчасти это может объясняться низкой временной локальностью. Однако в данном случае причина может также заключаться в том, что в данной реализации объем вычислений на одно обращение в память достаточно велик, что может приводить к недостаточной загруженности подсистемы памяти. В таком случае, несмотря на в целом неплохую эффективность работы с памятью, производительность будет достаточно низка.
2.3 Possible methods and considerations for parallel implementation of the algorithm
Итерация dqds практически полностью последовательна. Единственная возможность - одновременное выполнение операции умножения (2.3) и операции (2.4) умножения и сложения, что дает небольшой выигрыш в производительности.
2.4 Scalability of the algorithm and its implementations
Алгоритм не является масштабируемым, максимального эффекта ускорения можно добиться на двух независимых процессорах.
Проведём исследование масштабируемости вширь реализации согласно методике. Исследование проводилось на суперкомпьютере "Ломоносов-2" Суперкомпьютерного комплекса Московского университета.
Набор и границы значений изменяемых параметров запуска реализации алгоритма:
- число процессоров 1;
- размер матрицы [1000 : 260000] с шагом 500.
В результате проведённых экспериментов был получен следующий диапазон эффективности реализации алгоритма:
- минимальная эффективность реализации 2.559e-06%;
- максимальная эффективность реализации 2.492e-08%.
На следующих рисунках приведены графики производительности и эффективности выбранной реализации DQDS в зависимости от изменяемых параметров запуска.
Исследованная параллельная реализация на языке C
2.5 Dynamic characteristics and efficiency of the algorithm implementation
2.6 Conclusions for different classes of computer architecture
Эффективное выполнение алгоритма возможно только на вычислительных устройствах с одним или двумя ядрами.
2.7 Existing implementations of the algorithm
Сам алгоритм dqds реализован в функции xBDSQR пакета LAPACK и используется при её вызове без расчёта сингулярных векторов.
3 References
- ↑ 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.
