Уровень алгоритма

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

Материал из Алговики
Перейти к навигации Перейти к поиску
[непроверенная версия][досмотренная версия]
 
(не показано 35 промежуточных версий 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> (параметер итерации, см. [[#Математическое описание алгоритма|математическое описание алгоритма]]), отслеживание сходимости, а также применение различных оптимизационных "хитростей". Отметим, что внеитерационная часть алгоритма не существенна с точки зрения структуры вычислений, т.к. основные затраты ложаться на вычисления внутри итерации. Подробности и варианты внеитерационной части, а также анализ сходимости можно найти в соответствующей литературе
+
Вычислительным ядром алгоритма является именно внутренняя итерация, вне итераций происходит подбор сдвига <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>  
Строка 25: Строка 25:
 
=== Математическое описание алгоритма ===
 
=== Математическое описание алгоритма ===
  
 +
==== Вспомогательные сведения ====
 
Для понимания математических основ 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> - верхняя треугольная матрица с положительной диагональю.
Строка 39: Строка 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>q_j=a^2_j,\; e_j=b^2,</math> что в итоге приводит к так называемому алгоритму qds.
+
и вычислению <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.
 
Формулы алгоритма следующие:
 
Формулы алгоритма следующие:
  
Строка 54: Строка 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-итерации.
+
Такая математическая запись наиболее компактна и соответствует так называемой 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>
Строка 84: Строка 91:
 
=== Вычислительное ядро алгоритма ===
 
=== Вычислительное ядро алгоритма ===
  
Вычислительным ядром алгоритма является производимое n-1 раз вычисление, содержащие по одной операции сложения, вычитания и деления, а также две операции умножения.  
+
Вычислительным ядром алгоритма является последовательный расчёт квадратов диагональных (<math>\widehat{q}_j</math>) и внедиагональных (<math>\widehat{e}_k</math>) элементов выходной матрицы. Учитывая использование вспомогательных переменных расчёт каждой новой пары содержит по одной операции сложения, вычитания и деления, а также две операции умножения.
  
 
=== Макроструктура алгоритма ===
 
=== Макроструктура алгоритма ===
  
Алгоритм состоит из отдельного вычисления начального значения вспомогательной переменной <math>d</math> и последующим (n-1)-кратным выполнением повторяющейся последовательности из 5 операций (+,/,*,*,-).  
+
Алгоритм состоит из отдельного вычисления начального значения вспомогательной переменной <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>
Строка 149: Строка 158:
 
=== Особенности реализации последовательного алгоритма ===
 
=== Особенности реализации последовательного алгоритма ===
  
Алгоритм на языке Matlab может быть записать так:
+
Алгоритм на языке Matlab может быть записан так:
 
<source lang="matlab">
 
<source lang="matlab">
  
Строка 163: Строка 172:
 
</source>
 
</source>
  
 
+
Как говорилось в [[#Схема реализации последовательного алгоритма|cхеме реализации последовательного алгоритма]], вычисляемые данные записываются сразу на место входных.
=== Локальность данных и вычислений ===
 
 
 
Легко видеть, что локальность данных высока. Её легко повысить, размещая рядом соответствующие элементы массивов e и q, однако, это не оказывает существенного влияния на производительность в силу константного количества операций относительно объема обрабатываемых данных.
 
  
 
=== Возможные способы и особенности параллельной реализации алгоритма ===
 
=== Возможные способы и особенности параллельной реализации алгоритма ===
  
Итерация dqds практически полностью последовательна. Единственная возможность - одновременное выполнение операции умножения (2.3) и опрерации (2.4) умножения и сложения, что дает небольшой выигрыш в производительности.  
+
Итерация dqds практически полностью последовательна. Единственная возможность - одновременное выполнение операции умножения (2.3) и операции (2.4) умножения и сложения, что дает небольшой выигрыш в производительности.
 
 
=== Масштабируемость алгоритма и его реализации ===
 
  
Алгоритм не является масштабируемым, максимального эффекта ускорения можно добиться на двух независимых процессорах.
+
Сам алгоритм dqds реализован в функции xBDSQR пакета LAPACK и используется при её вызове без расчёта сингулярных векторов.
  
 +
=== Результаты прогонов ===
  
 
=== Выводы для классов архитектур ===
 
=== Выводы для классов архитектур ===
  
 
Эффективное выполнение алгоритма возможно только на вычислительных устройствах с одним или двумя ядрами.  
 
Эффективное выполнение алгоритма возможно только на вычислительных устройствах с одним или двумя ядрами.  
 
=== Существующие реализации алгоритма ===
 
Перечень реализаций приведен на странице [[Алгоритм dqds нахождения сингулярных чисел двухдиагональной матрицы#Существующие реализации алгоритма|алгоритма dqds]].
 
  
 
== Литература ==
 
== Литература ==
Строка 189: Строка 191:
  
 
[[Категория:Статьи в работе]]
 
[[Категория:Статьи в работе]]
 +
 +
[[En:One step of the dqds algorithm]]

Текущая версия на 09:49, 8 июля 2022



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


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

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] итерационно используя следующие три шага:

  1. Выбрать сдвиг [math]\tau_i[/math] меньший младшего собственного значения [math]T_i.[/math]
  2. Вычислить разложение Холецкого [math]T_i-\tau^2_iI=B_i^TB_i,[/math] где [math]B_i[/math] - верхняя треугольная матрица с положительной диагональю.
  3. [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 Вычислительное ядро алгоритма

Вычислительным ядром алгоритма является последовательный расчёт квадратов диагональных ([math]\widehat{q}_j[/math]) и внедиагональных ([math]\widehat{e}_k[/math]) элементов выходной матрицы. Учитывая использование вспомогательных переменных расчёт каждой новой пары содержит по одной операции сложения, вычитания и деления, а также две операции умножения.

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

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

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

Отметим, что выходные данные сразу могут быть записаны на место входных (это учтено в схеме), также для хранения вспомогательных переменных [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 Последовательная сложность алгоритма

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

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

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

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

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


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;

Как говорилось в cхеме реализации последовательного алгоритма, вычисляемые данные записываются сразу на место входных.

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

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

Сам алгоритм dqds реализован в функции xBDSQR пакета LAPACK и используется при её вызове без расчёта сингулярных векторов.

2.3 Результаты прогонов

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

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

3 Литература

  1. 1,0 1,1 1,2 Деммель Д. Вычислительная линейная алгебра. – М : Мир, 2001.
  2. Hogben L. (ed.). Handbook of linear algebra. – CRC Press, 2006.
  3. Fernando K. V., Parlett B. N. Accurate singular values and differential qd algorithms //Numerische Mathematik. – 1994. – Т. 67. – №. 2. – С. 191-229.
  4. Parlett B. N., Marques O. A. An implementation of the dqds algorithm (positive case) //Linear Algebra and its Applications. – 2000. – Т. 309. – №. 1. – С. 217-259.
  5. 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.