Итерация алгоритма dqds: различия между версиями
| [непроверенная версия] | [досмотренная версия] |
ASA (обсуждение | вклад) Нет описания правки |
|||
| (не показано 38 промежуточных версий 5 участников) | |||
| Строка 17: | Строка 17: | ||
'''Итерация алгоритма dqds''' является одним шагом алгоритма dqds нахождения сингулярных чисел двухдиагональной матрицы. | '''Итерация алгоритма dqds''' является одним шагом алгоритма dqds нахождения сингулярных чисел двухдиагональной матрицы. | ||
Сам алгоритм '''dqds''' (''differential quotient-difference algorithm with shifts'')<ref name="vla">Деммель Д. Вычислительная линейная алгебра. – М : Мир, 2001.</ref><ref name="hola">Hogben L. (ed.). Handbook of linear algebra. – CRC Press, 2006.</ref> строит последовательность | Сам алгоритм '''dqds''' (''differential quotient-difference algorithm with shifts'')<ref name="vla">Деммель Д. Вычислительная линейная алгебра. – М : Мир, 2001.</ref><ref name="hola">Hogben L. (ed.). Handbook of linear algebra. – CRC Press, 2006.</ref> строит последовательность двухдиагональных матриц, сходящуюся к диагональной матрице, содержащей квадраты искомых сингулярных чисел. Его особенностью является высокая точность решения задачи. | ||
Вычислительным ядром алгоритма является | Вычислительным ядром алгоритма является именно внутренняя итерация, вне итераций происходит подбор сдвига <math>\delta</math> (параметер итерации, см. [[#Математическое описание алгоритма|математическое описание алгоритма]]), отслеживание сходимости, а также применение различных оптимизационных "хитростей". Отметим, что внеитерационная часть алгоритма не существенна с точки зрения структуры вычислений, т.к. основные затраты ложатся на вычисления внутри итерации. Подробности и варианты внеитерационной части, а также анализ сходимости можно найти в соответствующей литературе | ||
<ref>Fernando K. V., Parlett B. N. Accurate singular values and differential qd algorithms //Numerische Mathematik. – 1994. – Т. 67. – №. 2. – С. 191-229.</ref> | <ref>Fernando K. V., Parlett B. N. Accurate singular values and differential qd algorithms //Numerische Mathematik. – 1994. – Т. 67. – №. 2. – С. 191-229.</ref> | ||
<ref>Parlett B. N., Marques O. A. An implementation of the dqds algorithm (positive case) //Linear Algebra and its Applications. – 2000. – Т. 309. – №. 1. – С. 217-259.</ref> | <ref>Parlett B. N., Marques O. A. An implementation of the dqds algorithm (positive case) //Linear Algebra and its Applications. – 2000. – Т. 309. – №. 1. – С. 217-259.</ref> | ||
<ref>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.</ref>. | <ref>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.</ref>. | ||
=== Математическое описание алгоритма === | === Математическое описание алгоритма === | ||
==== Вспомогательные сведения ==== | |||
Для понимания математических основ dqds-итерации полезно рассмотреть кратко её вывод, частично отражающий и историю возникновения алгоритма (подробности можно найти в <ref name="vla"/>). | Для понимания математических основ dqds-итерации полезно рассмотреть кратко её вывод, частично отражающий и историю возникновения алгоритма (подробности можно найти в <ref name="vla"/>). | ||
За основу dqds-алгоритма удобно взять так называемую LR-итерацию, предшествующую хорошо-известной QR-итерации. LR-алгоритм, начиная с входной матрицы <math>T_0>0,</math> строит сходящуюся последовательность подобных <math>T_0</math> матриц <math>T_i>0,</math> итерационно используя следующие три шага: | За основу 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>\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> - верхняя треугольная матрица с положительной диагональю. | ||
| Строка 58: | Строка 40: | ||
\widehat{b}^2_j = b^2_j\dot (a^2_{j+1}/a^2_j) | \widehat{b}^2_j = b^2_j\dot (a^2_{j+1}/a^2_j) | ||
</math> | </math> | ||
и вычислению <math>\widehat{a}^2_n = a^2_n-\widehat{b}^2_{n-1}-\delta.</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. | ||
Формулы алгоритма следующие: | Формулы алгоритма следующие: | ||
| Строка 73: | Строка 55: | ||
Здесь <math>q_j, \; j \in [1,n]</math> и <math>e_j, \; j \in [1,n-1]</math> - квадраты | Здесь <math>q_j, \; j \in [1,n]</math> и <math>e_j, \; j \in [1,n-1]</math> - квадраты элементов главной и верхней побочной диагонали соответственно. Крышка означает выходные переменные, а | ||
<math>\delta</math> - сдвиг (параметр алгоритма). | <math>\delta</math> - сдвиг (параметр алгоритма). | ||
Такая математическая запись наиболее компактна и | Такая математическая запись наиболее компактна и соответствует так называемой qds-итерации. | ||
==== Математическое описание итерации алгоритма dqds==== | |||
Представим теперь математическую запись, приближенную к dqds-итерации (с математической точки зрения qds и dqds-итерации эквивалентны) с введенными вспомогательными переменными | Представим теперь математическую запись, приближенную к dqds-итерации (с математической точки зрения qds и dqds-итерации эквивалентны) с введенными вспомогательными переменными | ||
<math>t_j</math> и <math>d_j:</math> | <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> | :<math> | ||
| Строка 103: | Строка 91: | ||
=== Вычислительное ядро алгоритма === | === Вычислительное ядро алгоритма === | ||
Вычислительным ядром алгоритма является | Вычислительным ядром алгоритма является последовательный расчёт квадратов диагональных (<math>\widehat{q}_j</math>) и внедиагональных (<math>\widehat{e}_k</math>) элементов выходной матрицы. Учитывая использование вспомогательных переменных расчёт каждой новой пары содержит по одной операции сложения, вычитания и деления, а также две операции умножения. | ||
=== Макроструктура алгоритма === | === Макроструктура алгоритма === | ||
Алгоритм состоит из отдельного вычисления начального значения вспомогательной переменной <math>d</math> | Алгоритм состоит из отдельного вычисления начального значения вспомогательной переменной <math>d,</math> последующего (n-1)-кратного выполнения повторяющейся последовательности из 5 операций (+,/,*,*,-) для вычисления квадратов диагональных (<math>\widehat{q}_j</math>) и внедиагональных (<math>\widehat{e}_k</math>) элементов выходной матрицы и завершающего вычисления крайнего значения <math>\widehat{q}_n</math>. | ||
=== Схема реализации последовательного алгоритма === | === Схема реализации последовательного алгоритма === | ||
Отметим, что выходные данные сразу могут быть записаны на место входных (это учтено в схеме), также для хранения вспомогательных переменных <math>t_j</math> и <math>d_j</math> достаточно двух перезаписываемых переменных. | Отметим, что выходные данные сразу могут быть записаны на место входных (это учтено в схеме), также для хранения вспомогательных переменных <math>t_j</math> и <math>d_j</math> достаточно двух перезаписываемых переменных. Таким образом элементы главной (<math>q_j</math>) и побочной (<math>e_k</math>) диагонали входной матрицы последовательно перезаписываются соответствующими элементами выходной матрицы. | ||
Последовательность исполнения метода следующая: | |||
1. Вычисляется начальное значение вспомогательной переменной <math>d = q_1-\delta.</math> | 1. Вычисляется начальное значение вспомогательной переменной <math>d = q_1-\delta.</math> | ||
| Строка 168: | Строка 158: | ||
=== Особенности реализации последовательного алгоритма === | === Особенности реализации последовательного алгоритма === | ||
Алгоритм на языке Matlab может быть | Алгоритм на языке Matlab может быть записан так: | ||
<source lang="matlab"> | <source lang="matlab"> | ||
| Строка 182: | Строка 172: | ||
</source> | </source> | ||
Как говорилось в [[#Схема реализации последовательного алгоритма|cхеме реализации последовательного алгоритма]], вычисляемые данные записываются сразу на место входных. | |||
=== Возможные способы и особенности параллельной реализации алгоритма === | === Возможные способы и особенности параллельной реализации алгоритма === | ||
Итерация dqds практически полностью последовательна. Единственная возможность - одновременное выполнение операции умножения (2.3) и | Итерация dqds практически полностью последовательна. Единственная возможность - одновременное выполнение операции умножения (2.3) и операции (2.4) умножения и сложения, что дает небольшой выигрыш в производительности. | ||
Сам алгоритм dqds реализован в функции xBDSQR пакета LAPACK и используется при её вызове без расчёта сингулярных векторов. | |||
=== Результаты прогонов === | |||
=== Выводы для классов архитектур === | === Выводы для классов архитектур === | ||
Эффективное выполнение алгоритма возможно только на вычислительных устройствах с одним или двумя ядрами. | Эффективное выполнение алгоритма возможно только на вычислительных устройствах с одним или двумя ядрами. | ||
== Литература == | == Литература == | ||
| Строка 208: | Строка 191: | ||
[[Категория:Статьи в работе]] | [[Категория:Статьи в работе]] | ||
[[En:One step of the dqds algorithm]] | |||
Текущая версия от 06:49, 8 июля 2022
| Алгоритм 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 нахождения сингулярных чисел двухдиагональной матрицы.
Сам алгоритм dqds (differential quotient-difference algorithm with shifts)[1][2] строит последовательность двухдиагональных матриц, сходящуюся к диагональной матрице, содержащей квадраты искомых сингулярных чисел. Его особенностью является высокая точность решения задачи. Вычислительным ядром алгоритма является именно внутренняя итерация, вне итераций происходит подбор сдвига [math]\displaystyle{ \delta }[/math] (параметер итерации, см. математическое описание алгоритма), отслеживание сходимости, а также применение различных оптимизационных "хитростей". Отметим, что внеитерационная часть алгоритма не существенна с точки зрения структуры вычислений, т.к. основные затраты ложатся на вычисления внутри итерации. Подробности и варианты внеитерационной части, а также анализ сходимости можно найти в соответствующей литературе [3] [4] [5].
Математическое описание алгоритма
Вспомогательные сведения
Для понимания математических основ dqds-итерации полезно рассмотреть кратко её вывод, частично отражающий и историю возникновения алгоритма (подробности можно найти в [1]). За основу dqds-алгоритма удобно взять так называемую LR-итерацию, предшествующую хорошо-известной QR-итерации. LR-алгоритм, начиная с входной симметричной и положительно определенной матрицы [math]\displaystyle{ T_0\gt 0, }[/math] строит сходящуюся последовательность подобных [math]\displaystyle{ T_0 }[/math] матриц [math]\displaystyle{ T_i\gt 0, }[/math] итерационно используя следующие три шага:
- Выбрать сдвиг [math]\displaystyle{ \tau_i }[/math] меньший младшего собственного значения [math]\displaystyle{ T_i. }[/math]
- Вычислить разложение Холецкого [math]\displaystyle{ T_i-\tau^2_iI=B_i^TB_i, }[/math] где [math]\displaystyle{ B_i }[/math] - верхняя треугольная матрица с положительной диагональю.
- [math]\displaystyle{ T_{i+1}=B_iB_i^T+\tau_i^2I. }[/math]
Отметим, что два шага LR-итерации с нулевым сдвигом эквивалентны одному шагу QR-итерации. Итерационная процедура приводит матрицу к диагональному виду, тем самым вычисляя собственные значения исходной матрицы. LR-алгоритм достаточно легко может быть переформулирован с заметными упрощениями для задачи поиска сингулярных значений двухдиагональных матриц. А именно, будем вычислять последовательность двухдиагональных матриц [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-1} }[/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-1}. }[/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{b}^2_{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{a}^2_n = a^2_n-\widehat{b}^2_{n-1}-\delta. }[/math] Очевидно, что работу с извлечением квадратных корней выгодно вести лишь после окончания работы алгоритма, поэтому можно ввести замену [math]\displaystyle{ q_j=a^2_j,\; e_j=b^2, }[/math] что в итоге приводит к так называемому алгоритму qds. Формулы алгоритма следующие:
- [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
Представим теперь математическую запись, приближенную к dqds-итерации (с математической точки зрения qds и dqds-итерации эквивалентны) с введенными вспомогательными переменными
[math]\displaystyle{ t_j }[/math] и [math]\displaystyle{ d_j. }[/math] Итерация алгоритма dqds преобразует входную двухдиагональную матрицу [math]\displaystyle{ B }[/math] в выходную [math]\displaystyle{ \widehat{B}. }[/math]
Входные и выходные данные: [math]\displaystyle{ q_{j}, \; j\in [1,n], \; e_{k}, \; k\in [1,n-1] }[/math] - квадраты элементов главной и побочной диагонали входной матрицы [math]\displaystyle{ B }[/math], [math]\displaystyle{ \widehat{q}_j , \; \widehat{e}_k }[/math] - то же для вычисляемой матрицы [math]\displaystyle{ \widehat{B}. }[/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{ \widehat{q}_j }[/math]) и внедиагональных ([math]\displaystyle{ \widehat{e}_k }[/math]) элементов выходной матрицы. Учитывая использование вспомогательных переменных расчёт каждой новой пары содержит по одной операции сложения, вычитания и деления, а также две операции умножения.
Макроструктура алгоритма
Алгоритм состоит из отдельного вычисления начального значения вспомогательной переменной [math]\displaystyle{ d, }[/math] последующего (n-1)-кратного выполнения повторяющейся последовательности из 5 операций (+,/,*,*,-) для вычисления квадратов диагональных ([math]\displaystyle{ \widehat{q}_j }[/math]) и внедиагональных ([math]\displaystyle{ \widehat{e}_k }[/math]) элементов выходной матрицы и завершающего вычисления крайнего значения [math]\displaystyle{ \widehat{q}_n }[/math].
Схема реализации последовательного алгоритма
Отметим, что выходные данные сразу могут быть записаны на место входных (это учтено в схеме), также для хранения вспомогательных переменных [math]\displaystyle{ t_j }[/math] и [math]\displaystyle{ d_j }[/math] достаточно двух перезаписываемых переменных. Таким образом элементы главной ([math]\displaystyle{ q_j }[/math]) и побочной ([math]\displaystyle{ e_k }[/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 реализация вычисления позволяет достичь высокой точности[1].
Последовательная сложность алгоритма
Для выполнения одной итерации dqds необходимо выполнить:
- [math]\displaystyle{ n-1 }[/math] делений,
- [math]\displaystyle{ 2n-2 }[/math] умножений,
- [math]\displaystyle{ 2n-1 }[/math] сложений/вычитаний.
Таким образом одна dqds-итерация имеет линейную сложность.
Информационный граф
Ресурс параллелизма алгоритма
Как видно из информационного графа алгоритма, на каждом шаге основного цикла возможно лишь параллельное выполнение операции умножения (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;
Как говорилось в cхеме реализации последовательного алгоритма, вычисляемые данные записываются сразу на место входных.
Возможные способы и особенности параллельной реализации алгоритма
Итерация dqds практически полностью последовательна. Единственная возможность - одновременное выполнение операции умножения (2.3) и операции (2.4) умножения и сложения, что дает небольшой выигрыш в производительности.
Сам алгоритм dqds реализован в функции xBDSQR пакета LAPACK и используется при её вызове без расчёта сингулярных векторов.
Результаты прогонов
Выводы для классов архитектур
Эффективное выполнение алгоритма возможно только на вычислительных устройствах с одним или двумя ядрами.
Литература
- ↑ 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.