Участник:Davletshuna Alexandra/Метод Якоби вычисления сингулярных чисел и векторов: различия между версиями
Frolov (обсуждение | вклад) |
|||
(не показано 50 промежуточных версий 3 участников) | |||
Строка 1: | Строка 1: | ||
+ | {{Assignment|ASA|Frolov}} | ||
+ | |||
{{algorithm | {{algorithm | ||
Строка 9: | Строка 11: | ||
}} | }} | ||
− | Авторы статьи: [[Участница:Davletshuna Alexandra|Давлетшина Александра]] (группа 619) | + | Авторы статьи: [[Участница:Davletshuna Alexandra|Давлетшина Александра]] (группа 619); [[Участница:Alek sandra|Зайцева Александра]]. Вклад каждого участника считать равноценным. |
== Свойства и структура алгоритма == | == Свойства и структура алгоритма == | ||
Строка 38: | Строка 40: | ||
симметричную матрицу <math>A = (G^TG)</math> можно привести к диагональному виду: | симметричную матрицу <math>A = (G^TG)</math> можно привести к диагональному виду: | ||
<math>RG^TGR^T</math>, где | <math>RG^TGR^T</math>, где | ||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
<math> | <math> | ||
Строка 75: | Строка 71: | ||
Обозначим <math>s = \sin \theta</math> и <math>c = \cos \theta</math>. | Обозначим <math>s = \sin \theta</math> и <math>c = \cos \theta</math>. | ||
− | Перемножив матрицы <math>RG^TGR^T</math>, получим следующие выражения: | + | Перемножив матрицы <math>RG^TGR^T</math>, получим следующие выражения<ref>http://www.astro.tsu.ru/OsChMet/5_2.html</ref>: |
<math>\begin{align} | <math>\begin{align} | ||
Строка 92: | Строка 88: | ||
<math> | <math> | ||
\begin{align} | \begin{align} | ||
− | \tau &= {a_{jj} - a_{kk}}/{2a_{jk}} \\ | + | \tau &= ({a_{jj} - a_{kk}})/{2a_{jk}} \\ |
t &= {sign(\tau)}/(|\tau| + \sqrt{1 + \tau^2}) \\ | t &= {sign(\tau)}/(|\tau| + \sqrt{1 + \tau^2}) \\ | ||
c &= 1/(\sqrt{1 + t^2}) \\ | c &= 1/(\sqrt{1 + t^2}) \\ | ||
Строка 99: | Строка 95: | ||
Остановка вычислений алгоритма наступает, когда внедиагональные элементы становятся меньше некоторого наперед заданного числа <math>\varepsilon</math>. | Остановка вычислений алгоритма наступает, когда внедиагональные элементы становятся меньше некоторого наперед заданного числа <math>\varepsilon</math>. | ||
+ | |||
+ | В случае, когда <math> a_{jj} = a_{kk} </math> считаем <math> sign(\tau) = 1 </math>, т. е. <math> \theta = | ||
+ | \pi/4 </math>. | ||
=== Вычислительное ядро алгоритма === | === Вычислительное ядро алгоритма === | ||
− | Вычислительное ядро алгоритма составляет непосредственно вычисление одностороннего вращения Якоби, т.е. вычисление элементов матрицы <math>g_{jm}^{(i+1)} = g_{mj}^{(i+1)}</math> и <math>g_{km}^{(i+1)} = g_{mk}^{(i+1)}</math>, <math>m \ne j,k</math>, а также вычисления элементов <math>g_{jj}^{(i+1)} </math> и <math>g_{kk}^{(i+1)} </math>, получившихся после умножения промежуточной матрицы на матрицу вращений Якоби | + | Вычислительное ядро алгоритма составляет непосредственно вычисление одностороннего вращения Якоби, т.е. вычисление элементов матрицы <math>g_{jm}^{(i+1)} = g_{mj}^{(i+1)}</math> и <math>g_{km}^{(i+1)} = g_{mk}^{(i+1)}</math>, <math>m \ne j,k</math>, а также вычисления элементов <math>g_{jj}^{(i+1)} </math> и <math>g_{kk}^{(i+1)} </math>, получившихся после умножения промежуточной матрицы на матрицу вращений Якоби: |
+ | |||
+ | <math>\begin{align} | ||
+ | g_{jj}^{(i+1)} &= c^2\, g_{jj}^{(i)} - 2\, s c \,g_{jk}^{(i)} + s^2\, g_{kk}^{(i)} \\ | ||
+ | g_{kk}^{(i+1)} &= s^2 \,g_{jj}^{(i)} + 2 s c\, g_{jk}^{(i)} + c^2 \, g_{kk}^{(i)} \\ | ||
+ | g_{jk}^{(i+1)} &= g_{kj}^{(i+1)} = (c^2 - s^2 ) \, g_{jk}^{(i)} + s c \, (g_{kk}^{(i)} - g_{jj}^{(i)} ) \\ | ||
+ | g_{jm}^{(i+1)} &= g_{mj}^{(i+1)} = c \, g_{jm}^{(i)} - s \, g_{km}^{(i)} & m \ne j,k \\ | ||
+ | g_{km}^{(i+1)} &= g_{mk}^{(i+1)} = s \, g_{jm}^{(i)} + c \, g_{km}^{(i)} & m \ne j,k \\ | ||
+ | g_{ml}^{(i+1)} &= g_{ml}^{(i)} &m,l \ne j,k | ||
+ | \end{align}</math> | ||
=== Макроструктура алгоритма === | === Макроструктура алгоритма === | ||
Строка 139: | Строка 147: | ||
<math> | <math> | ||
\begin{align} | \begin{align} | ||
− | \tau &= {a_{jj} - a_{kk}}/{2a_{jk}} \\ | + | \tau &= ({a_{jj} - a_{kk}})/{2a_{jk}} \\ |
t &= {sign(\tau)}/(|\tau| + \sqrt{1 + \tau^2}) \\ | t &= {sign(\tau)}/(|\tau| + \sqrt{1 + \tau^2}) \\ | ||
c &= 1/(\sqrt{1 + t^2}) \\ | c &= 1/(\sqrt{1 + t^2}) \\ | ||
Строка 150: | Строка 158: | ||
=== Последовательная сложность алгоритма === | === Последовательная сложность алгоритма === | ||
+ | Для одного вращения Якоби необходимо выполнить <math>O(n)</math> операций умножения и сложения. | ||
+ | |||
+ | Всего вариантов выбора <math>j</math> и <math>k</math> существует <math>\frac{n(n-1)}{2}</math>. | ||
+ | |||
+ | Таким образом последовательный алгоритм Якоби вычисления сингулярных чисел и векторов является алгоритмом с кубической сложностью. | ||
=== Информационный граф === | === Информационный граф === | ||
+ | Граф алгоритма состоит из четырех групп вершин | ||
+ | |||
+ | (А)Первая группа вершин расположена в двумерной области, соответствующая ей операция - вычисление элементов <math>a_{jj},a_{jk},a_{kk}</math> матрицы произведения <math>A = (G^TG)</math>. Где <math>j</math> изменяется от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения, <math>k</math> изменяется от <math>j+1</math> до <math>n</math>, принимая все целочисленные значения. | ||
+ | |||
+ | Граф вычисления каждого из элементов приведен на рисунке 1, где значениям <math>g_{ij}</math> соответствуют обозначения <math>x_{ij}</math>, а элементам <math>g^T_{ij}</math> обозначения <math>z_{ij}</math>. Значение <math>s_{ij}</math> является значением результирующего элемента. | ||
+ | |||
+ | (Б)Вторая группа вершин расположена в двумерной области, соответствующая ей операция – вычисление <math>c</math> и <math>s</math>. | ||
+ | |||
+ | Информационный граф с входными и выходными данными вычисления пункта (Б) представлен на рисунке 2. | ||
+ | |||
+ | (В)Третья группа вершин расположена в двумерной области, соответствующая ей операция вычисление значений элементов <math>g_{jm}^{(i+1)} = g_{mj}^{(i+1)}</math> и <math>g_{km}^{(i+1)} = g_{mk}^{(i+1)}</math>, <math>m \ne j,k</math>. | ||
+ | |||
+ | (Г)Четвертая группа вершин расположена в двумерной области, соответствующая ей операция вычисление значений элементов <math>g_{jj}^{(i+1)} </math> и <math>g_{kk}^{(i+1)} </math>. | ||
+ | |||
+ | Элементы групп (В) и (Г) вычисляются аналогично вычислению элементов матрицы произведения, т.к. фактически представляют собой перемножение матриц <math>G</math> и <math>R</math> | ||
+ | |||
+ | [[File:Jacobi_2.jpg|none|thumb|1000px|Рисунок 1. Информационный граф с входными и выходными параметрами вычисления каждого элемента матрицы произведения.]] | ||
+ | |||
+ | [[File:Jacobi_1.jpg|none|thumb|1000px|Рисунок 2. Информационный граф с входными и выходными параметрами вычисления значений <math>c</math> и <math>s</math>.]] | ||
+ | |||
+ | Обобщенный граф метода Якоби представлен на рисунке 3, где буквами А, Б, В и Г обозначено вычисление групп вершин (А), (Б), (В) и (Г). | ||
+ | |||
+ | [[File:Graf21.png|none|thumb|2000px|Рисунок 3. Обобщенный граф алгоритма вычисления сингулярных чисел и векторов методом Якоби.]] | ||
=== Ресурс параллелизма алгоритма === | === Ресурс параллелизма алгоритма === | ||
Строка 158: | Строка 194: | ||
# вычислить <math>t</math> | # вычислить <math>t</math> | ||
# вычислить <math>c</math> | # вычислить <math>c</math> | ||
+ | # вычислить <math>s</math> | ||
# выполнить поворот матрицы <math>G</math> с параллельным выполнением <math>2(n-2)</math> операций вычисления элементов матрицы <math>g_{jm}^{(i+1)} = g_{mj}^{(i+1)}</math> и <math>g_{km}^{(i+1)} = g_{mk}^{(i+1)}</math>, <math>m \ne j,k</math>, а также вычисления элементов <math>g_{jj}^{(i+1)} </math> и <math>g_{kk}^{(i+1)} </math>. | # выполнить поворот матрицы <math>G</math> с параллельным выполнением <math>2(n-2)</math> операций вычисления элементов матрицы <math>g_{jm}^{(i+1)} = g_{mj}^{(i+1)}</math> и <math>g_{km}^{(i+1)} = g_{mk}^{(i+1)}</math>, <math>m \ne j,k</math>, а также вычисления элементов <math>g_{jj}^{(i+1)} </math> и <math>g_{kk}^{(i+1)} </math>. | ||
Строка 191: | Строка 228: | ||
=== Масштабируемость алгоритма и его реализации === | === Масштабируемость алгоритма и его реализации === | ||
+ | Исследование масштабируемости параллельной реализации метода Якоби вычисления сингулярных чисел и векторов проводилось на тестовой средe суперкомпьютера "Ломоносов", с использованием компилятора ICC и библиотеки impi/4.0.3. Алгоритм реализован на языке программирования C с использованием механизмов MPI. | ||
+ | |||
+ | Алгоритм работает по следующей схеме: | ||
+ | #На вход всем процессорам подается матрица <math> G </math>. | ||
+ | #Для всех <math> j </math> и <math> k </math> таких, что <math> j \in [1,n-1] </math>, <math> k \in [j+1,n] </math> | ||
+ | ##Вычисляется элемент <math> a_{jk} </math> путем перемножения матриц <math> G^TG </math>, где на каждом процессоре вычисляется своя часть матрицы. После чего результат отправляется на все процессоры. | ||
+ | ##Если полученный элемент оказывается больше некоторого значения, то начинается последовательная часть программы: Вычисление элементов <math> \tau, t, c, s </math>. | ||
+ | ##Затем, выполняется операция поворота, путем перемножения матриц <math> G </math> и <math> R </math>. Данная операция выполняется параллельно на нескольких процессорах, где каждый процессор вычисляет свою часть матрицы. После чего результат пересылается всем процессорам. | ||
+ | |||
+ | Данная схема повторяется до тех пор, пока матрица <math> A = G^TG </math> не станет достаточно близкой к диагональной. | ||
+ | |||
+ | Для проверки этого факта, в конце каждой итерации цикла вычисляется матрица <math> A </math>. Данная операция выполняется параллельно на нескольких процессорах, где каждый процессор вычисляет свою часть матрицы, после чего результат пересылается всем процессорам. | ||
+ | |||
+ | |||
+ | Набор и границы значений изменяемых параметров запуска реализации алгоритма: | ||
+ | |||
+ | * число процессоров [1 : 128] с шагом равным степени двойки; | ||
+ | * размер матрицы [256 : 1792] с шагом 256. | ||
+ | |||
+ | В результате проведённых экспериментов был получен следующий диапазон эффективности реализации алгоритма: | ||
+ | |||
+ | * минимальная эффективность реализации 0,01%; | ||
+ | * максимальная эффективность реализации 16,52%. | ||
+ | |||
+ | На следующих рисунках приведены графики производительности и эффективности выбранной реализации метода Якоби в зависимости от изменяемых параметров запуска. | ||
+ | |||
+ | [[File:Performance.jpg|none|thumb|1000px|Рисунок 1. Параллельная реализация метода Якоби. Изменение производительности в зависимости от числа процессоров и размера матрицы.]] | ||
+ | |||
+ | [[File:Эффективность.png|none|thumb|1000px|Рисунок 1. Параллельная реализация метода Якоби. Изменение эффективности в зависимости от числа процессоров и размера матрицы.]] | ||
+ | |||
+ | Построим оценки масштабируемости данной реализации метода Якоби: | ||
+ | * По числу процессов: при увеличении числа процессов эффективность на рассмотренной области изменений параметров запуска сначала увеличивается (самый пик достигается при 16 процессах), а затем начинает снижаться. Это характеризуется тем, что накладные расходы на параллельное выполнение (общение между процессами) начинает превалировать над вычислениями. | ||
+ | |||
+ | *По размеру задачи: при увеличении размера задачи эффективность возрастает. Это подтверждает предположение о том, что размер задачи сильно влияет на эффективность выполнения приложения. Оценка показывает, что с ростом размера задачи эффективность на рассмотренной области значений параметров запуска сильно увеличивается. | ||
+ | |||
+ | *По двум направлениям: при рассмотрении увеличения как вычислительной сложности, так и числа процессов на всей рассмотренной области значений эффективность увеличивается, но медленно и с перепадами. | ||
+ | |||
+ | [https://www.dropbox.com/home?preview=mainmpi.c Исследованная параллельная реализация метода Якоби] | ||
=== Динамические характеристики и эффективность реализации алгоритма === | === Динамические характеристики и эффективность реализации алгоритма === | ||
Строка 197: | Строка 272: | ||
=== Существующие реализации алгоритма === | === Существующие реализации алгоритма === | ||
− | Метод Якоби нахождения сингулярных значений и векторов реализован в таких пакетах, как LAPACK, LINPACK/EISPACK, Intel MKL <ref> | + | Метод Якоби нахождения сингулярных значений и векторов реализован в таких пакетах, как LAPACK <ref>http://www.netlib.org/scalapack/slug/</ref>, LINPACK/EISPACK <ref> http://www.netlib.org/linpack/ </ref>, Intel MKL <ref>https://software.intel.com/ru-ru/node/521153#CAF8CA08-1C69-48F8-B0C3-B2B27AFA92FA</ref>. |
Так же функция Якоби реализована в пакете Matlab. | Так же функция Якоби реализована в пакете Matlab. |
Текущая версия на 17:04, 17 декабря 2016
Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Frolov и ASA. |
Метод Якоби вычисления сингулярных чисел и векторов | |
Последовательный алгоритм | |
Последовательная сложность | [math]O(n^3)[/math] |
Объём входных данных | [math]n^2[/math] |
Объём выходных данных | [math]2n^2 + n[/math] |
Параллельный алгоритм | |
Высота ярусно-параллельной формы | [math]O(n^2)[/math] |
Ширина ярусно-параллельной формы | [math]O(n)[/math] |
Авторы статьи: Давлетшина Александра (группа 619); Зайцева Александра. Вклад каждого участника считать равноценным.
Содержание
- 1 Свойства и структура алгоритма
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 2.1 Особенности реализации последовательного алгоритма
- 2.2 Локальность данных и вычислений
- 2.3 Возможные способы и особенности параллельной реализации алгоритма
- 2.4 Масштабируемость алгоритма и его реализации
- 2.5 Динамические характеристики и эффективность реализации алгоритма
- 2.6 Выводы для классов архитектур
- 2.7 Существующие реализации алгоритма
- 3 Литература
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
Сингулярным разложением называется разложение прямоугольной вещественной или комплексной матрицы, имеющее широкое применение, в силу своей наглядной геометрической интерпретации, при решении многих прикладных задач. Сингулярное разложение матрицы А позволяет вычислять сингулярные числа данной матрицы, а также, левые и правые сингулярные векторы матрицы А:
• левые сингулярные векторы матрицы А — это собственные векторы матрицы [math]AA^*[/math];
• правые сингулярные векторы матрицы A — это собственные векторы матрицы [math]A^*A[/math].
Геометрическая интерпретация сингулярного разложения заключается в следующем факте из геометрии: Образом любого линейного преобразования, заданного с помощью матрицы [math]m[/math]x[math]n[/math], примененного к единичной сфере является гиперэллипсоид.
Из многочисленных разложений матриц, наиболее удобным является сингулярное разложение матрицы в виде: [math]A=U\Sigma V^*[/math], где [math]U, V[/math] – унитарные матрицы, а [math]\Sigma[/math] – диагональная матрица с вещественными положительными числами на диагонали [1].
Диагональные элементы матрицы [math]\Sigma[/math] называются сингулярными числами матрицы [math]A[/math], а столбцы матриц [math]U,V[/math] левыми и правыми сингулярными векторами соответственно.
Алгоритм Якоби сингулярного разложения матрицы был предложен одним из первых в 1846 году. Он приводит прямоугольную матрицу к диагональной матрице с помощью последовательности элементарных вращений. Метод может найти все сингулярные значения с очень высокой точностью. Однако его производительность является довольно низкой, в сравнении с конкурирующими методами.
Неявный[2] метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице [math]A = G^TG[/math]. Это значит, что на каждом шаге вычисляется вращение Якоби [math]J[/math], с помощью которого матрица [math]G^TG[/math] неявно пересчитывется в [math]J^TG^TGJ[/math]; вращение выбрано так, чтобы пара внедиагональных элементов из [math]G^TG[/math] обратилась в нули в матрице [math]J^TG^TGJ[/math]. При этом ни [math]G^TG[/math], ни [math]J^TG^TGJ[/math] не вычисляются в явном виде; вместо них вычисляется матрица [math]GJ[/math]. Поэтому алгоритм называется методом односторонних вращений.
1.2 Математическое описание алгоритма
Исходные данные: матрица [math]G[/math] (элементы [math]g_{ij}, i, j = 1, \ldots, n[/math]).
Вычисляемые данные: матрица [math]\Sigma = diag(\sigma_{i})[/math], где [math]\sigma_{i}[/math] - сингулярные числа, [math]U[/math] - матрица левых сингулярных векторов, [math]V[/math] - матрица правых сингулярных векторов.
симметричную матрицу [math]A = (G^TG)[/math] можно привести к диагональному виду: [math]RG^TGR^T[/math], где
[math] R(j,k,\theta) = \begin{matrix} \\ \\ \\ j \\ \\ k \\ \\ \\ \\ \end{matrix} [/math] [math] \begin{bmatrix} & 1 & & & & & & & & \\ & & 1 & & & & & & & \\ & & & \ddots & & & & & & \\ & & & & \cos(\theta) & & -\sin(\theta) & & & \\ & & & & & \ddots & & & & \\ & & & & \sin(\theta) & & \cos(\theta) & & & \\ & & & & & & & \ddots & & \\ & & & & & & & & 1 & \\ & & & & & & & & & 1 \\ \end{bmatrix} [/math]
Обозначим [math]s = \sin \theta[/math] и [math]c = \cos \theta[/math].
Перемножив матрицы [math]RG^TGR^T[/math], получим следующие выражения[3]:
[math]\begin{align} a_{jj}^{(i+1)} &= c^2\, a_{jj}^{(i)} - 2\, s c \,a_{jk}^{(i)} + s^2\, a_{kk}^{(i)} \\ a_{kk}^{(i+1)} &= s^2 \,a_{jj}^{(i)} + 2 s c\, a_{jk}^{(i)} + c^2 \, a_{kk}^{(i)} \\ a_{jk}^{(i+1)} &= a_{kj}^{(i+1)} = (c^2 - s^2 ) \, a_{jk}^{(i)} + s c \, (a_{kk}^{(i)} - a_{jj}^{(i)} ) \\ a_{jm}^{(i+1)} &= a_{mj}^{(i+1)} = c \, a_{jm}^{(i)} - s \, a_{km}^{(i)} & m \ne j,k \\ a_{km}^{(i+1)} &= a_{mk}^{(i+1)} = s \, a_{jm}^{(i)} + c \, a_{km}^{(i)} & m \ne j,k \\ a_{ml}^{(i+1)} &= a_{ml}^{(i)} &m,l \ne j,k \end{align}[/math]
Т.к внедиагональные элементы должны быть равны 0, то выбираем [math]\theta[/math] так, чтобы [math]a_{jk}^{(i+1)} = 0[/math] и [math]a_{kj}^{(i+1)} = 0[/math].
Тогда решая уравнение получаем, что
[math] \begin{align} \tau &= ({a_{jj} - a_{kk}})/{2a_{jk}} \\ t &= {sign(\tau)}/(|\tau| + \sqrt{1 + \tau^2}) \\ c &= 1/(\sqrt{1 + t^2}) \\ s &= ct \\ \end{align}[/math]
Остановка вычислений алгоритма наступает, когда внедиагональные элементы становятся меньше некоторого наперед заданного числа [math]\varepsilon[/math].
В случае, когда [math] a_{jj} = a_{kk} [/math] считаем [math] sign(\tau) = 1 [/math], т. е. [math] \theta = \pi/4 [/math].
1.3 Вычислительное ядро алгоритма
Вычислительное ядро алгоритма составляет непосредственно вычисление одностороннего вращения Якоби, т.е. вычисление элементов матрицы [math]g_{jm}^{(i+1)} = g_{mj}^{(i+1)}[/math] и [math]g_{km}^{(i+1)} = g_{mk}^{(i+1)}[/math], [math]m \ne j,k[/math], а также вычисления элементов [math]g_{jj}^{(i+1)} [/math] и [math]g_{kk}^{(i+1)} [/math], получившихся после умножения промежуточной матрицы на матрицу вращений Якоби:
[math]\begin{align} g_{jj}^{(i+1)} &= c^2\, g_{jj}^{(i)} - 2\, s c \,g_{jk}^{(i)} + s^2\, g_{kk}^{(i)} \\ g_{kk}^{(i+1)} &= s^2 \,g_{jj}^{(i)} + 2 s c\, g_{jk}^{(i)} + c^2 \, g_{kk}^{(i)} \\ g_{jk}^{(i+1)} &= g_{kj}^{(i+1)} = (c^2 - s^2 ) \, g_{jk}^{(i)} + s c \, (g_{kk}^{(i)} - g_{jj}^{(i)} ) \\ g_{jm}^{(i+1)} &= g_{mj}^{(i+1)} = c \, g_{jm}^{(i)} - s \, g_{km}^{(i)} & m \ne j,k \\ g_{km}^{(i+1)} &= g_{mk}^{(i+1)} = s \, g_{jm}^{(i)} + c \, g_{km}^{(i)} & m \ne j,k \\ g_{ml}^{(i+1)} &= g_{ml}^{(i)} &m,l \ne j,k \end{align}[/math]
1.4 Макроструктура алгоритма
Основную часть метода составляет процедура односторонних вращений матрицы, которая в свою очередь состоит из:
- Вычисление [math]\tau, c, s, t[/math].
- Одностороннее вращение матрицы.
Что требует [math]O(n)[/math] операций.
После этого необходимо:
- Вычислить сингулярные числа [math]\sigma_i[/math]. Что требует [math]O(n)[/math] операций.
- Вычислить [math]U[/math] и [math]V[/math]. Что требует [math]O(n^2)[/math] операций.
1.5 Схема реализации последовательного алгоритма
Применение метода Якоби к матрице [math]A = (G^TG)[/math] можно описать следующим образом[4]:
Если [math]G = U\Sigma V^T[/math] и [math]\Sigma = diag(\sigma_{i})[/math], то алгоритм вычисляет сингулярные числа [math]\sigma_{i}[/math], матрицу [math]U[/math] левых сингулярных векторов и матрицу [math]V[/math] правых сингулярных векторов по следующей схеме:
repeat for [math]j=1[/math] to [math]n-1[/math] for [math]k=j+1[/math] to [math]n[/math] обратиться к процедуре One-Sided-Jacobi-Rotation[math](G,j,k)[/math] end for end for пока [math]G^TG[/math] не станет достаточно близка к диагональной матрице Положить [math]\sigma_{i} = ||G(:,i)||_2 [/math] (норма [math]i[/math]-го столбца в [math]G[/math]) Положить [math]U = [u_{1},\dots, u_{n}][/math], где [math]u_{i} = G(:,i)/\sigma_{i}[/math] Положить [math]V = J[/math], где [math]J[/math] - накопленное произведение вращений Якоби.
Где процедура One-Sided-Jacobi-Rotation является метод односторонних вращений Якоби:
proc One-Sided-Jacobi-Rotation[math](G,j,k)[/math] вычислить [math]a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk}[/math] и [math]a_{kk} = (G^TG)_{kk}[/math] if [math]|a_{jk}|[/math] не слишком мал [math] \begin{align} \tau &= ({a_{jj} - a_{kk}})/{2a_{jk}} \\ t &= {sign(\tau)}/(|\tau| + \sqrt{1 + \tau^2}) \\ c &= 1/(\sqrt{1 + t^2}) \\ s &= ct \\ G &= GR(j,k,\theta) \ \dots c=cos\theta, s=sin\theta \end{align}[/math] if нужны правые сингулярные векторы [math]J = JR(j,k,\theta)[/math] end if end if
1.6 Последовательная сложность алгоритма
Для одного вращения Якоби необходимо выполнить [math]O(n)[/math] операций умножения и сложения.
Всего вариантов выбора [math]j[/math] и [math]k[/math] существует [math]\frac{n(n-1)}{2}[/math].
Таким образом последовательный алгоритм Якоби вычисления сингулярных чисел и векторов является алгоритмом с кубической сложностью.
1.7 Информационный граф
Граф алгоритма состоит из четырех групп вершин
(А)Первая группа вершин расположена в двумерной области, соответствующая ей операция - вычисление элементов [math]a_{jj},a_{jk},a_{kk}[/math] матрицы произведения [math]A = (G^TG)[/math]. Где [math]j[/math] изменяется от [math]1[/math] до [math]n-1[/math], принимая все целочисленные значения, [math]k[/math] изменяется от [math]j+1[/math] до [math]n[/math], принимая все целочисленные значения.
Граф вычисления каждого из элементов приведен на рисунке 1, где значениям [math]g_{ij}[/math] соответствуют обозначения [math]x_{ij}[/math], а элементам [math]g^T_{ij}[/math] обозначения [math]z_{ij}[/math]. Значение [math]s_{ij}[/math] является значением результирующего элемента.
(Б)Вторая группа вершин расположена в двумерной области, соответствующая ей операция – вычисление [math]c[/math] и [math]s[/math].
Информационный граф с входными и выходными данными вычисления пункта (Б) представлен на рисунке 2.
(В)Третья группа вершин расположена в двумерной области, соответствующая ей операция вычисление значений элементов [math]g_{jm}^{(i+1)} = g_{mj}^{(i+1)}[/math] и [math]g_{km}^{(i+1)} = g_{mk}^{(i+1)}[/math], [math]m \ne j,k[/math].
(Г)Четвертая группа вершин расположена в двумерной области, соответствующая ей операция вычисление значений элементов [math]g_{jj}^{(i+1)} [/math] и [math]g_{kk}^{(i+1)} [/math].
Элементы групп (В) и (Г) вычисляются аналогично вычислению элементов матрицы произведения, т.к. фактически представляют собой перемножение матриц [math]G[/math] и [math]R[/math]
Обобщенный граф метода Якоби представлен на рисунке 3, где буквами А, Б, В и Г обозначено вычисление групп вершин (А), (Б), (В) и (Г).
1.8 Ресурс параллелизма алгоритма
Для выполнения одной итерации метода Якоби требуется последовательно выполнить несколько ярусов:
- вычислить [math]\tau[/math]
- вычислить [math]t[/math]
- вычислить [math]c[/math]
- вычислить [math]s[/math]
- выполнить поворот матрицы [math]G[/math] с параллельным выполнением [math]2(n-2)[/math] операций вычисления элементов матрицы [math]g_{jm}^{(i+1)} = g_{mj}^{(i+1)}[/math] и [math]g_{km}^{(i+1)} = g_{mk}^{(i+1)}[/math], [math]m \ne j,k[/math], а также вычисления элементов [math]g_{jj}^{(i+1)} [/math] и [math]g_{kk}^{(i+1)} [/math].
Суммарное количество итераций равно [math]\frac{n(n-1)}{2}[/math] (Различные варианты выбора [math]j[/math] и [math]k[/math]).
Асимптотически метод сходится квадратично.
Таким образом, при классификации по высоте ЯПФ метод Якоби относится к алгоритмам со сложностью [math]O(n^2)[/math] , а при классификации по ширине ЯПФ — к алгоритмам со сложностью [math]O(n)[/math].
1.9 Входные и выходные данные алгоритма
Входные данные: плотная матрица [math]G[/math] размера [math]n[/math]х[math]n[/math], где [math]A= G^TG[/math] – симметрическая матрица.
Объём входных данных: [math]n^2[/math] (т.к. для хранения матрицы [math]G[/math] используем двумерный массив размера [math]n[/math]х[math]n[/math])
Выходные данные: [math]\sigma_{i}[/math] - сингулярные числа, матрица [math]U[/math] левых сингулярных векторов и матрица [math]V[/math] правых сингулярных векторов.
Объём выходных данных: [math]2n^2 + n[/math] (т.к. необходимо хранить вектор сингулярных чисел длины [math]n[/math], а также две матрицы [math]U, V[/math] левых и правых сингулярных векторов размера [math]n[/math]х[math]n[/math]).
1.10 Свойства алгоритма
Метод Якоби является самым медленным из имеющихся алгоритмов вычисления сингулярных чисел и сингулярных векторов матрицы. Тем не менее, для матриц [math]G[/math], которые могут быть представлены в виде [math]G=DX[/math], где [math]D[/math] - диагональная, а [math]X[/math] - хорошо обусловлена, он способен вычислять сингулярные числа и вектора намного точнее, чем другие методы[5].
Вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных является линейной.
Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности (число итераций зависит от входных данных и порогового значения).
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
2.2 Локальность данных и вычислений
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
Исследование масштабируемости параллельной реализации метода Якоби вычисления сингулярных чисел и векторов проводилось на тестовой средe суперкомпьютера "Ломоносов", с использованием компилятора ICC и библиотеки impi/4.0.3. Алгоритм реализован на языке программирования C с использованием механизмов MPI.
Алгоритм работает по следующей схеме:
- На вход всем процессорам подается матрица [math] G [/math].
- Для всех [math] j [/math] и [math] k [/math] таких, что [math] j \in [1,n-1] [/math], [math] k \in [j+1,n] [/math]
- Вычисляется элемент [math] a_{jk} [/math] путем перемножения матриц [math] G^TG [/math], где на каждом процессоре вычисляется своя часть матрицы. После чего результат отправляется на все процессоры.
- Если полученный элемент оказывается больше некоторого значения, то начинается последовательная часть программы: Вычисление элементов [math] \tau, t, c, s [/math].
- Затем, выполняется операция поворота, путем перемножения матриц [math] G [/math] и [math] R [/math]. Данная операция выполняется параллельно на нескольких процессорах, где каждый процессор вычисляет свою часть матрицы. После чего результат пересылается всем процессорам.
Данная схема повторяется до тех пор, пока матрица [math] A = G^TG [/math] не станет достаточно близкой к диагональной.
Для проверки этого факта, в конце каждой итерации цикла вычисляется матрица [math] A [/math]. Данная операция выполняется параллельно на нескольких процессорах, где каждый процессор вычисляет свою часть матрицы, после чего результат пересылается всем процессорам.
Набор и границы значений изменяемых параметров запуска реализации алгоритма:
- число процессоров [1 : 128] с шагом равным степени двойки;
- размер матрицы [256 : 1792] с шагом 256.
В результате проведённых экспериментов был получен следующий диапазон эффективности реализации алгоритма:
- минимальная эффективность реализации 0,01%;
- максимальная эффективность реализации 16,52%.
На следующих рисунках приведены графики производительности и эффективности выбранной реализации метода Якоби в зависимости от изменяемых параметров запуска.
Построим оценки масштабируемости данной реализации метода Якоби:
- По числу процессов: при увеличении числа процессов эффективность на рассмотренной области изменений параметров запуска сначала увеличивается (самый пик достигается при 16 процессах), а затем начинает снижаться. Это характеризуется тем, что накладные расходы на параллельное выполнение (общение между процессами) начинает превалировать над вычислениями.
- По размеру задачи: при увеличении размера задачи эффективность возрастает. Это подтверждает предположение о том, что размер задачи сильно влияет на эффективность выполнения приложения. Оценка показывает, что с ростом размера задачи эффективность на рассмотренной области значений параметров запуска сильно увеличивается.
- По двум направлениям: при рассмотрении увеличения как вычислительной сложности, так и числа процессов на всей рассмотренной области значений эффективность увеличивается, но медленно и с перепадами.
Исследованная параллельная реализация метода Якоби
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
Метод Якоби нахождения сингулярных значений и векторов реализован в таких пакетах, как LAPACK [6], LINPACK/EISPACK [7], Intel MKL [8].
Так же функция Якоби реализована в пакете Matlab.
Т. к. алгоритм считается медленным, он не включен во многие известные пакеты.
3 Литература
1. Дж. Деммель Вычислительная линейная алгебра. Изд. Мир, 2001.
- ↑ http://algorithms.parallel.ru/
- ↑ Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)
- ↑ http://www.astro.tsu.ru/OsChMet/5_2.html
- ↑ Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-263)
- ↑ Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)
- ↑ http://www.netlib.org/scalapack/slug/
- ↑ http://www.netlib.org/linpack/
- ↑ https://software.intel.com/ru-ru/node/521153#CAF8CA08-1C69-48F8-B0C3-B2B27AFA92FA