Алговики - Вклад [ru]

Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов

2016-11-19T12:05:27Z

Светлана: /* Масштабируемость реализации алгоритма */

{{algorithm
| name = Метод Якоби вычисления сингулярных чисел и векторов
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math>O(n^2)</math>
| pf_width = <math>O(n)</math>
| input_data = <math>n^2</math>
| output_data = <math>2n^2 + n</math>
}}
Авторы статьи: [[Участник:Почернина Елена|Почернина Елена]], группа 601 (разделы [[#Общее описание алгоритма|1.1]], [[#Математическое описание алгоритма|1.2]], [[#Макроструктура алгоритма|1.4]], [[#Схема реализации последовательного алгоритма|1.5]], [[#Последовательная сложность алгоритма|1.6]], [[#Ресурс параллелизма алгоритма|1.8]], [[#Масштабируемость алгоритма и его реализации|2.4]], [[#Существующие реализации алгоритма|2.7]]), [[Участница:Светлана|Костюкова Светлана]], группа 601 (разделы [[#Вычислительное ядро алгоритма|1.3]], [[#Схема реализации последовательного алгоритма|1.5]], [[#Информационный граф|1.7]], [[#Входные и выходные данные алгоритма|1.9]], [[#Свойства алгоритма|1.10]], [[#Масштабируемость алгоритма и его реализации|2.4]])

== Свойства и структура алгоритма ==

=== Общее описание алгоритма ===

'''Метод Якоби''' позволяет найти сингулярное разложение плотной матрицы <math>G</math> размера <math>n</math>×<math>n</math>. Такое разложение можно записать в виде <math>G = U\Sigma V^*</math>, где <math>U, V</math> - унитарные (''в вещественном случае - ортогональные'') матрицы размера <math>n</math>×<math>n</math>, <math>\Sigma</math> - диагональная матрица размера <math>n</math>×<math>n</math> с вещественными положительными числами на главной диагонали. Диагональные элементы матрицы <math>\Sigma</math> называются сингулярными числами матрицы <math>G</math>, а столбцы матриц <math>U,V</math> - левыми и правыми сингулярными векторами.

''Неявный''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref> метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице <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>. Поэтому алгоритм называется ''методом односторонних вращений''.

Метод Якоби для вычисления сингулярных чисел и векторов является самым медленным из имеющихся алгоритмов поиска сингулярных чисел и векторов, но тем не менее, интерес к нему сохраняется. Для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. В частности, метод Якоби вычисляет сингулярные числа матрицы <math>G</math> с высокой точностью, если G может быть представлена в виде <math>G = DX</math>, где <math>D</math> – диагональная матрица, а <math>X</math> – хорошо обусловлена. В этом случае заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, в то время как другие алгоритмы включают в себя приведение матрицы к двухдиагональной форме, из-за чего и теряют все верные разряды во всех сингулярных числах, кроме старшего.

=== Математическое описание алгоритма ===
Исходные данные: матрица <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> - матрица правых сингулярных векторов.

Основной алгоритм заключается в следующем<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref>:

Для всех элементов матрицы <math>A^{(i)} = (G^TG)^{(i)}</math>, находящихся вне главной диагонали <math>(a_{12}, a_{13},\dots, a_{1n}, a_{23}, a_{24},\dots, a_{2n},\dots, a_{n-1n})</math> проверяем условие <math>|A_{jk}|<\varepsilon</math>:
Если <math>|A_{jk}^{(i)}| \le \varepsilon</math>, то переходим к следующему элементу,
Если <math>|A_{jk}^{(i)}| > \varepsilon</math>, то нужно выполнить вращение в плоскости <math>(j,k)</math>, после чего перейти к следующей итерации алгоритма.
Если на очередной итерации не было совершено ни одного вращения, значит матрица <math>A</math> стала достаточна близка к диагональной, и можно перейти к вычислению искомых значений:

<math>\sigma_i = ||G(:,i)||_2</math>;

<math>U = [u_1,\dots,u_n]</math>, где <math>u_i = G(:,i)/\sigma_i</math>;

<math>V = J</math>, где <math>J</math> - накопленное произведение вращений Якоби.

Одностороннее вращение в плоскости <math>(j,k)</math> осуществляет процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''':

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 &= \frac{a_{jj} - a_{kk}}{2a_{jk}} \\
t &= \frac{sign(\tau)}{|\tau| + \sqrt{1 + \tau^2}} \\
c &= \frac{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

Здесь <math>R(j,k,\theta)</math> - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 244-245)</ref>.

                                                                                
<math>
\begin{matrix}
& & & & j & & k & & & \\
\end{matrix}
</math>

<math>
J_i = 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>a_{jk}</math> не слишком мал"''. Подходящим критерием останова является условие
<math>\begin{align}
|a_{jk}| \ge \varepsilon\sqrt{a_{jj}a_{kk}}
\end{align}</math>.

=== Вычислительное ядро алгоритма ===

Вычислительное ядро составляют множественные вычисления элементов матрицы <math>a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk}</math> и <math>a_{kk} = (G^TG)_{kk}</math>, а также вычисление промежуточной матрицы <math>GR(j, k, \theta)</math> путем умножения исходной матрицы на матрицу вращения Якоби в процессе выполнения процедуры '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''.

После приведения матрицы <math>G^TG</math> к виду, достаточно близкому к диагональному, вычисляются сингулярные числа:

<math>\sigma_i = ||G(:,i)||_2</math>

Матрица левых сингулярных векторов <math>U = [u_1,..., u_n]</math>:

<math>u_i = G(:,i)/\sigma_i</math>

И (если требуется) матрица правых сингулярных векторов <math>V</math> путем накопления произведения вращений Якоби:

<math>
V = J,
J = JR(j, k, \theta)
</math>

=== Макроструктура алгоритма ===
Основную часть метода составляет процедура односторонних вращений матрицы (процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''). Рассмотрим подробнее структуру этой процедуры.

- Вычисление <math>a_{jk}, a_{jj}, a_{kk} (a_{ij} = (G^TG)_{ij})</math>. В связи с тем, что нет необходимости каждый раз вычислять произведение <math>G^TG</math>, для вычисления каждого элемента требуется <math>n</math> операций умножения и <math>(n-1)</math> операций сложения.

- Вычисление <math>\tau, c, s, t</math> не представляют значительной вычислительной сложности.

- Одностороннее вращение матрицы (умножение входной матрицы на матрицу вращения) требует порядка <math>O(n)</math> операций.

В итоге, процедура '''''One-Sided-Jacobi-Rotation''''' требует <math>O(n)</math> операций.

После выполнения алгоритма односторонних вращений необходимо вычислить искомые матрицы <math>S, U, V</math> (матрица сингулярных чисел и матрицы левых и правых векторов).

- Вычисление сингулярных чисел <math>\sigma_i</math> требует <math>n</math> операций умножения.

- Вычисление <math>U</math> требует <math>O(n^2)</math> операций. Аналогично для матрицы <math>V</math>.

=== Схема реализации последовательного алгоритма ===
Метод можно описать следующим образом<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-263)</ref>:

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> - накопленное произведение вращений Якоби.

Метод Якоби для нахождения сингулярных чисел и векторов матрицы <math>G</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 &= \frac{a_{jj} - a_{kk}}{2a_{jk}} \\
t &= \frac{sign(\tau)}{|\tau| + \sqrt{1 + \tau^2}} \\
c &= \frac{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

=== Последовательная сложность алгоритма ===
Для вычисления сингулярных чисел и векторов матрицы порядка <math>n</math> в последовательном варианте на каждой итерации алгоритма вызывается процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''', сложность которой равна <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>j</math> — меняется в диапазоне от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения;
* <math>k</math> — меняется в диапазоне от <math>j+1</math> до <math>n</math>, принимая все целочисленные значения.

'''Вторая''' группа вершин соответствует вычислению значений <math>c</math> и <math>s</math> для матрицы вращения

'''Третья''' группа вершин соответствует вычислению произведений матриц <math>GR(k,j, \theta)</math> и <math>JR(k,j, \theta)</math>

<center>
<div class="thumb">
<div class="thumbinner" style="width:{{#expr: 3 * (400 + 35) + 4 * (3 - 1) + 8}}px">
<gallery widths=400px heights=400px>
File:JacobiSVG.png|Рисунок 1. Граф одной итерации алгоритма без отображения входных и выходных данных.
File:Akj.png|Рисунок 2. Граф вычисления значений элементов <math>a_{jj}, a_{jk}, a_{kk}</math>.
</gallery>
<gallery widths=400px heights=1200px>
File:JacobiSVD2.png|Рисунок 3. Граф одной итерации алгоритма с вычислением матрицы вращения <math>R</math> отображением входных и выходных данных.
</gallery>
<gallery widths=400px heights=400px>
File:Gij.png|Рисунок 4. Граф вычисления значений элементов матрицы <math>GR</math>
File:Jij.png|Рисунок 5. Граф вычисления значений элементов матрицы <math>JR</math>
</gallery>
</div>
</div>
</center>

=== Ресурс параллелизма алгоритма ===
В данном алгоритме существует возможность параллельного вычисления вращений Якоби<ref>Halil Snopce, Ilir Spahiu Parallel Computation of the SVD, pp. 458-460, http://www.wseas.us/e-library/conferences/2011/Paris/ECC/ECC-73.pdf</ref><ref>V. Strumpen, H. Hoffmann, A. Agarwal A Stream Algorithm for the SVD, pp. 6-7, 23-24, http://publications.csail.mit.edu/lcs/pubs/pdf/MIT-LCS-TM-641.pdf</ref>. Максимальное количество вращений в данном методе равно <math>\frac{n(n-1)}{2}</math>. В параллельной реализации, если <math>n</math> - четное, можно выполнить <math>\frac{n}{2}</math> вращений параллельно, если <math>n</math> - нечетное, то <math>\frac{n-1}{2}</math> вращений.

Основную сложность в процедуре вращения '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''' представляет вычисление элементов <math>a_{jj}, a_{jk}, a_{kk}</math>, которое можно производить параллельно. Вычисление элементов <math>a_{jj}, a_{jk}, a_{kk}</math> требует порядка <math>O(logN)</math> операций. Вычисление элементов матриц <math>GR, JR</math> в параллельной реализации требует <math>O(1)</math> операций.

Таким образом, при классификации по высоте ЯПФ метод Якоби относится к алгоритмам со сложностью <math>O(n^2)</math>, а при классификации по ширине ЯПФ — к алгоритмам со сложностью <math>O(n)</math>.

=== Входные и выходные данные алгоритма ===
'''Входные данные''': плотная матрица <math>G</math> (элементы <math>g_{ij}</math>).
Дополнительные ограничения:
* <math>A= G^TG</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

'''Объём входных данных''': <math>n^2</math>

'''Выходные данные''': матрица <math>\Sigma = diag(\sigma_{i})</math>, где <math>\sigma_{i}</math> - сингулярные числа, матрица <math>U</math> левых сингулярных векторов и матрица <math>V</math> правых сингулярных векторов.

'''Объём выходных данных''': <math>2n^2 + n</math> (так как матрица сингулярных чисел <math>\Sigma</math> - диагональная, то достаточно хранить вектор сингулярных чисел (элементы <math>\sigma_i</math>), а также две матрицы <math>U, V</math> левых и правых сингулярных векторов).

=== Свойства алгоритма ===
Метод Якоби является самым медленным из имеющихся алгоритмов вычисления сингулярных чисел и сингулярных векторов матрицы. Тем не менее, для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. Кроме того, в методе Якоби заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, за счет этого метод достигает высокой относительной точности.

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов для метода Якоби является линейным, вычислительная мощность алгоритма также линейна.

Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности (число итераций зависит от входных данных и порогового значения).

== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===

=== Локальность данных и вычислений ===

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

=== Масштабируемость алгоритма и его реализации ===

==== Масштабируемость алгоритма ====
Оценки предела масштабируемости алгоритма были рассмотрены в разделе [[#Ресурс параллелизма алгоритма|Ресурс параллелизма алгоритма]].

==== Масштабируемость реализации алгоритма ====
Проведём исследование масштабируемости параллельной реализации метода Якоби нахождения сингулярных чисел и векторов. Исследование проводилось на суперкомпьютере "Ломоносов" [http://parallel.ru/cluster Суперкомпьютерного комплекса Московского университета].

Набор и границы значений изменяемых реализации алгоритма:

* число процессоров [1:16] с шагом 1;
* размер матрицы [100 : 900] с шагом 100.

При запуске использовалась версия mkl 11.2.0.

Компилятор icc запускался со следующими опциями: -DMKL_ILP64 -qopenmp -I${MKLROOT}/include -Wl,--start-group ${MKLROOT}/lib/intel64/libmkl_intel_ilp64.a ${MKLROOT}/lib/intel64/libmkl_core.a ${MKLROOT}/lib/intel64/libmkl_intel_thread.a -Wl,--end-group -lpthread -lm -ldl

На следующем рисунке приведен график среднего времени работы алгоритма в зависимости от изменяемых параметров запуска.

[[File:Averagetime.PNG|center|thumb|900px|Рисунок 6. Параллельная реализация метода Якоби вычисления сингулярных чисел и векторов. Изменение времени выполнения в зависимости от числа процессоров и размера матрицы]]
[[File:JacobiPerformance.PNG|center|thumb|900px|Рисунок 7. Параллельная реализация метода Якоби вычисления сингулярных чисел и векторов. Изменение производительности в зависимости от числа процессов и размера матрицы]]

В результате проведенных экспериментов был получен следующий диапазон [[Глоссарий#Эффективность реализации|эффективности реализации]] алгоритма:

* минимальная эффективность реализации: 0.0095%
* максимальная эффективность реализации: 0.929%

При увеличении числа процессов эффективность на рассмотренной области изменений параметров запуска падает.

[https://software.intel.com/ru-ru/node/521152#248AD4D8-005A-4800-9A73-428952BC1E44 Исследованная параллельная реализация, реализованная в библиотеке Intel MKL]

=== Динамические характеристики и эффективность реализации алгоритма ===

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

=== Существующие реализации алгоритма ===
Несмотря на то, что метод Якоби является самым медленным из методов вычисления сингулярных чисел и векторов матрицы, он реализован в некоторых пакетах, например, LAPACK, ScaLAPACK, Intel MKL<ref>Intel® Math Kernel Library, https://software.intel.com/ru-ru/node/521153#CAF8CA08-1C69-48F8-B0C3-B2B27AFA92FA</ref>.

== Литература ==

Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов

2016-11-19T11:41:55Z

Светлана: /* Масштабируемость реализации алгоритма */

{{algorithm
| name = Метод Якоби вычисления сингулярных чисел и векторов
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math>O(n^2)</math>
| pf_width = <math>O(n)</math>
| input_data = <math>n^2</math>
| output_data = <math>2n^2 + n</math>
}}
Авторы статьи: [[Участник:Почернина Елена|Почернина Елена]], группа 601 (разделы [[#Общее описание алгоритма|1.1]], [[#Математическое описание алгоритма|1.2]], [[#Макроструктура алгоритма|1.4]], [[#Схема реализации последовательного алгоритма|1.5]], [[#Последовательная сложность алгоритма|1.6]], [[#Ресурс параллелизма алгоритма|1.8]], [[#Масштабируемость алгоритма и его реализации|2.4]], [[#Существующие реализации алгоритма|2.7]]), [[Участница:Светлана|Костюкова Светлана]], группа 601 (разделы [[#Вычислительное ядро алгоритма|1.3]], [[#Схема реализации последовательного алгоритма|1.5]], [[#Информационный граф|1.7]], [[#Входные и выходные данные алгоритма|1.9]], [[#Свойства алгоритма|1.10]], [[#Масштабируемость алгоритма и его реализации|2.4]])

== Свойства и структура алгоритма ==

=== Общее описание алгоритма ===

'''Метод Якоби''' позволяет найти сингулярное разложение плотной матрицы <math>G</math> размера <math>n</math>×<math>n</math>. Такое разложение можно записать в виде <math>G = U\Sigma V^*</math>, где <math>U, V</math> - унитарные (''в вещественном случае - ортогональные'') матрицы размера <math>n</math>×<math>n</math>, <math>\Sigma</math> - диагональная матрица размера <math>n</math>×<math>n</math> с вещественными положительными числами на главной диагонали. Диагональные элементы матрицы <math>\Sigma</math> называются сингулярными числами матрицы <math>G</math>, а столбцы матриц <math>U,V</math> - левыми и правыми сингулярными векторами.

''Неявный''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref> метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице <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>. Поэтому алгоритм называется ''методом односторонних вращений''.

Метод Якоби для вычисления сингулярных чисел и векторов является самым медленным из имеющихся алгоритмов поиска сингулярных чисел и векторов, но тем не менее, интерес к нему сохраняется. Для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. В частности, метод Якоби вычисляет сингулярные числа матрицы <math>G</math> с высокой точностью, если G может быть представлена в виде <math>G = DX</math>, где <math>D</math> – диагональная матрица, а <math>X</math> – хорошо обусловлена. В этом случае заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, в то время как другие алгоритмы включают в себя приведение матрицы к двухдиагональной форме, из-за чего и теряют все верные разряды во всех сингулярных числах, кроме старшего.

=== Математическое описание алгоритма ===
Исходные данные: матрица <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> - матрица правых сингулярных векторов.

Основной алгоритм заключается в следующем<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref>:

Для всех элементов матрицы <math>A^{(i)} = (G^TG)^{(i)}</math>, находящихся вне главной диагонали <math>(a_{12}, a_{13},\dots, a_{1n}, a_{23}, a_{24},\dots, a_{2n},\dots, a_{n-1n})</math> проверяем условие <math>|A_{jk}|<\varepsilon</math>:
Если <math>|A_{jk}^{(i)}| \le \varepsilon</math>, то переходим к следующему элементу,
Если <math>|A_{jk}^{(i)}| > \varepsilon</math>, то нужно выполнить вращение в плоскости <math>(j,k)</math>, после чего перейти к следующей итерации алгоритма.
Если на очередной итерации не было совершено ни одного вращения, значит матрица <math>A</math> стала достаточна близка к диагональной, и можно перейти к вычислению искомых значений:

<math>\sigma_i = ||G(:,i)||_2</math>;

<math>U = [u_1,\dots,u_n]</math>, где <math>u_i = G(:,i)/\sigma_i</math>;

<math>V = J</math>, где <math>J</math> - накопленное произведение вращений Якоби.

Одностороннее вращение в плоскости <math>(j,k)</math> осуществляет процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''':

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 &= \frac{a_{jj} - a_{kk}}{2a_{jk}} \\
t &= \frac{sign(\tau)}{|\tau| + \sqrt{1 + \tau^2}} \\
c &= \frac{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

Здесь <math>R(j,k,\theta)</math> - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 244-245)</ref>.

                                                                                
<math>
\begin{matrix}
& & & & j & & k & & & \\
\end{matrix}
</math>

<math>
J_i = 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>a_{jk}</math> не слишком мал"''. Подходящим критерием останова является условие
<math>\begin{align}
|a_{jk}| \ge \varepsilon\sqrt{a_{jj}a_{kk}}
\end{align}</math>.

=== Вычислительное ядро алгоритма ===

Вычислительное ядро составляют множественные вычисления элементов матрицы <math>a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk}</math> и <math>a_{kk} = (G^TG)_{kk}</math>, а также вычисление промежуточной матрицы <math>GR(j, k, \theta)</math> путем умножения исходной матрицы на матрицу вращения Якоби в процессе выполнения процедуры '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''.

После приведения матрицы <math>G^TG</math> к виду, достаточно близкому к диагональному, вычисляются сингулярные числа:

<math>\sigma_i = ||G(:,i)||_2</math>

Матрица левых сингулярных векторов <math>U = [u_1,..., u_n]</math>:

<math>u_i = G(:,i)/\sigma_i</math>

И (если требуется) матрица правых сингулярных векторов <math>V</math> путем накопления произведения вращений Якоби:

<math>
V = J,
J = JR(j, k, \theta)
</math>

=== Макроструктура алгоритма ===
Основную часть метода составляет процедура односторонних вращений матрицы (процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''). Рассмотрим подробнее структуру этой процедуры.

- Вычисление <math>a_{jk}, a_{jj}, a_{kk} (a_{ij} = (G^TG)_{ij})</math>. В связи с тем, что нет необходимости каждый раз вычислять произведение <math>G^TG</math>, для вычисления каждого элемента требуется <math>n</math> операций умножения и <math>(n-1)</math> операций сложения.

- Вычисление <math>\tau, c, s, t</math> не представляют значительной вычислительной сложности.

- Одностороннее вращение матрицы (умножение входной матрицы на матрицу вращения) требует порядка <math>O(n)</math> операций.

В итоге, процедура '''''One-Sided-Jacobi-Rotation''''' требует <math>O(n)</math> операций.

После выполнения алгоритма односторонних вращений необходимо вычислить искомые матрицы <math>S, U, V</math> (матрица сингулярных чисел и матрицы левых и правых векторов).

- Вычисление сингулярных чисел <math>\sigma_i</math> требует <math>n</math> операций умножения.

- Вычисление <math>U</math> требует <math>O(n^2)</math> операций. Аналогично для матрицы <math>V</math>.

=== Схема реализации последовательного алгоритма ===
Метод можно описать следующим образом<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-263)</ref>:

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> - накопленное произведение вращений Якоби.

Метод Якоби для нахождения сингулярных чисел и векторов матрицы <math>G</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 &= \frac{a_{jj} - a_{kk}}{2a_{jk}} \\
t &= \frac{sign(\tau)}{|\tau| + \sqrt{1 + \tau^2}} \\
c &= \frac{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

=== Последовательная сложность алгоритма ===
Для вычисления сингулярных чисел и векторов матрицы порядка <math>n</math> в последовательном варианте на каждой итерации алгоритма вызывается процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''', сложность которой равна <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>j</math> — меняется в диапазоне от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения;
* <math>k</math> — меняется в диапазоне от <math>j+1</math> до <math>n</math>, принимая все целочисленные значения.

'''Вторая''' группа вершин соответствует вычислению значений <math>c</math> и <math>s</math> для матрицы вращения

'''Третья''' группа вершин соответствует вычислению произведений матриц <math>GR(k,j, \theta)</math> и <math>JR(k,j, \theta)</math>

<center>
<div class="thumb">
<div class="thumbinner" style="width:{{#expr: 3 * (400 + 35) + 4 * (3 - 1) + 8}}px">
<gallery widths=400px heights=400px>
File:JacobiSVG.png|Рисунок 1. Граф одной итерации алгоритма без отображения входных и выходных данных.
File:Akj.png|Рисунок 2. Граф вычисления значений элементов <math>a_{jj}, a_{jk}, a_{kk}</math>.
</gallery>
<gallery widths=400px heights=1200px>
File:JacobiSVD2.png|Рисунок 3. Граф одной итерации алгоритма с вычислением матрицы вращения <math>R</math> отображением входных и выходных данных.
</gallery>
<gallery widths=400px heights=400px>
File:Gij.png|Рисунок 4. Граф вычисления значений элементов матрицы <math>GR</math>
File:Jij.png|Рисунок 5. Граф вычисления значений элементов матрицы <math>JR</math>
</gallery>
</div>
</div>
</center>

=== Ресурс параллелизма алгоритма ===
В данном алгоритме существует возможность параллельного вычисления вращений Якоби<ref>Halil Snopce, Ilir Spahiu Parallel Computation of the SVD, pp. 458-460, http://www.wseas.us/e-library/conferences/2011/Paris/ECC/ECC-73.pdf</ref><ref>V. Strumpen, H. Hoffmann, A. Agarwal A Stream Algorithm for the SVD, pp. 6-7, 23-24, http://publications.csail.mit.edu/lcs/pubs/pdf/MIT-LCS-TM-641.pdf</ref>. Максимальное количество вращений в данном методе равно <math>\frac{n(n-1)}{2}</math>. В параллельной реализации, если <math>n</math> - четное, можно выполнить <math>\frac{n}{2}</math> вращений параллельно, если <math>n</math> - нечетное, то <math>\frac{n-1}{2}</math> вращений.

Основную сложность в процедуре вращения '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''' представляет вычисление элементов <math>a_{jj}, a_{jk}, a_{kk}</math>, которое можно производить параллельно. Вычисление элементов <math>a_{jj}, a_{jk}, a_{kk}</math> требует порядка <math>O(logN)</math> операций. Вычисление элементов матриц <math>GR, JR</math> в параллельной реализации требует <math>O(1)</math> операций.

Таким образом, при классификации по высоте ЯПФ метод Якоби относится к алгоритмам со сложностью <math>O(n^2)</math>, а при классификации по ширине ЯПФ — к алгоритмам со сложностью <math>O(n)</math>.

=== Входные и выходные данные алгоритма ===
'''Входные данные''': плотная матрица <math>G</math> (элементы <math>g_{ij}</math>).
Дополнительные ограничения:
* <math>A= G^TG</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

'''Объём входных данных''': <math>n^2</math>

'''Выходные данные''': матрица <math>\Sigma = diag(\sigma_{i})</math>, где <math>\sigma_{i}</math> - сингулярные числа, матрица <math>U</math> левых сингулярных векторов и матрица <math>V</math> правых сингулярных векторов.

'''Объём выходных данных''': <math>2n^2 + n</math> (так как матрица сингулярных чисел <math>\Sigma</math> - диагональная, то достаточно хранить вектор сингулярных чисел (элементы <math>\sigma_i</math>), а также две матрицы <math>U, V</math> левых и правых сингулярных векторов).

=== Свойства алгоритма ===
Метод Якоби является самым медленным из имеющихся алгоритмов вычисления сингулярных чисел и сингулярных векторов матрицы. Тем не менее, для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. Кроме того, в методе Якоби заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, за счет этого метод достигает высокой относительной точности.

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов для метода Якоби является линейным, вычислительная мощность алгоритма также линейна.

Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности (число итераций зависит от входных данных и порогового значения).

== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===

=== Локальность данных и вычислений ===

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

=== Масштабируемость алгоритма и его реализации ===

==== Масштабируемость алгоритма ====
Оценки предела масштабируемости алгоритма были рассмотрены в разделе [[#Ресурс параллелизма алгоритма|Ресурс параллелизма алгоритма]].

==== Масштабируемость реализации алгоритма ====
Проведём исследование масштабируемости параллельной реализации метода Якоби нахождения сингулярных чисел и векторов. Исследование проводилось на суперкомпьютере "Ломоносов" [http://parallel.ru/cluster Суперкомпьютерного комплекса Московского университета].

Набор и границы значений изменяемых реализации алгоритма:

* число процессоров [1:16] с шагом 1;
* размер матрицы [100 : 900] с шагом 100.

При запуске использовалась версия mkl 11.2.0.

Компилятор icc запускался со следующими опциями: -DMKL_ILP64 -qopenmp -I${MKLROOT}/include -Wl,--start-group ${MKLROOT}/lib/intel64/libmkl_intel_ilp64.a ${MKLROOT}/lib/intel64/libmkl_core.a ${MKLROOT}/lib/intel64/libmkl_intel_thread.a -Wl,--end-group -lpthread -lm -ldl

На следующем рисунке приведен график среднего времени работы алгоритма в зависимости от изменяемых параметров запуска.

[[File:Averagetime.PNG|center|thumb|900px|Рисунок 6. Параллельная реализация метода Якоби вычисления сингулярных чисел и векторов. Изменение времени выполнения в зависимости от числа процессоров и размера матрицы]]
[[File:JacobiPerformance.PNG|center|thumb|900px|Рисунок 7. Параллельная реализация метода Якоби вычисления сингулярных чисел и векторов. Изменение производительности в зависимости от числа процессов и размера матрицы]]

При увеличении числа процессов эффективность на рассмотренной области изменений параметров запуска падает.

[https://software.intel.com/ru-ru/node/521152#248AD4D8-005A-4800-9A73-428952BC1E44 Исследованная параллельная реализация, реализованная в библиотеке Intel MKL]

=== Динамические характеристики и эффективность реализации алгоритма ===

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

=== Существующие реализации алгоритма ===
Несмотря на то, что метод Якоби является самым медленным из методов вычисления сингулярных чисел и векторов матрицы, он реализован в некоторых пакетах, например, LAPACK, ScaLAPACK, Intel MKL<ref>Intel® Math Kernel Library, https://software.intel.com/ru-ru/node/521153#CAF8CA08-1C69-48F8-B0C3-B2B27AFA92FA</ref>.

== Литература ==

Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов

2016-11-15T20:22:06Z

Светлана: /* Масштабируемость реализации алгоритма */

{{Assignment|ASA}}
{{Assignment}}
{{algorithm
| name = Метод Якоби вычисления сингулярных чисел и векторов
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math>O(n^2)</math>
| pf_width = <math>O(n)</math>
| input_data = <math>n^2</math>
| output_data = <math>2n^2 + n</math>
}}
Авторы статьи: [[Участник:Почернина Елена|Почернина Елена]], группа 601 (разделы [[#Общее описание алгоритма|1.1]], [[#Математическое описание алгоритма|1.2]], [[#Макроструктура алгоритма|1.4]], [[#Схема реализации последовательного алгоритма|1.5]], [[#Последовательная сложность алгоритма|1.6]], [[#Ресурс параллелизма алгоритма|1.8]], [[#Масштабируемость алгоритма и его реализации|2.4]], [[#Существующие реализации алгоритма|2.7]]), [[Участница:Светлана|Костюкова Светлана]], группа 601 (разделы [[#Вычислительное ядро алгоритма|1.3]], [[#Схема реализации последовательного алгоритма|1.5]], [[#Информационный граф|1.7]], [[#Входные и выходные данные алгоритма|1.9]], [[#Свойства алгоритма|1.10]], [[#Масштабируемость алгоритма и его реализации|2.4]])

== Свойства и структура алгоритма ==

=== Общее описание алгоритма ===

'''Метод Якоби''' позволяет найти сингулярное разложение плотной матрицы <math>G</math> размера <math>n</math>×<math>n</math>. Такое разложение можно записать в виде <math>G = U\Sigma V^*</math>, где <math>U, V</math> - унитарные (''в вещественном случае - ортогональные'') матрицы размера <math>n</math>×<math>n</math>, <math>\Sigma</math> - диагональная матрица размера <math>n</math>×<math>n</math> с вещественными положительными числами на главной диагонали. Диагональные элементы матрицы <math>\Sigma</math> называются сингулярными числами матрицы <math>G</math>, а столбцы матриц <math>U,V</math> - левыми и правыми сингулярными векторами.

''Неявный''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref> метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице <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>. Поэтому алгоритм называется ''методом односторонних вращений''.

Метод Якоби для вычисления сингулярных чисел и векторов является самым медленным из имеющихся алгоритмов поиска сингулярных чисел и векторов, но тем не менее, интерес к нему сохраняется. Для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. В частности, метод Якоби вычисляет сингулярные числа матрицы <math>G</math> с высокой точностью, если G может быть представлена в виде <math>G = DX</math>, где <math>D</math> – диагональная матрица, а <math>X</math> – хорошо обусловлена. В этом случае заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, в то время как другие алгоритмы включают в себя приведение матрицы к двухдиагональной форме, из-за чего и теряют все верные разряды во всех сингулярных числах, кроме старшего.

=== Математическое описание алгоритма ===
Исходные данные: матрица <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> - матрица правых сингулярных векторов.

Основной алгоритм заключается в следующем<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref>:

Для всех элементов матрицы <math>A^{(i)} = (G^TG)^{(i)}</math>, находящихся вне главной диагонали <math>(a_{12}, a_{13},\dots, a_{1n}, a_{23}, a_{24},\dots, a_{2n},\dots, a_{n-1n})</math> проверяем условие <math>|A_{jk}|<\varepsilon</math>:
Если <math>|A_{jk}^{(i)}| \le \varepsilon</math>, то переходим к следующему элементу,
Если <math>|A_{jk}^{(i)}| > \varepsilon</math>, то нужно выполнить вращение в плоскости <math>(j,k)</math>, после чего перейти к следующей итерации алгоритма.
Если на очередной итерации не было совершено ни одного вращения, значит матрица <math>A</math> стала достаточна близка к диагональной, и можно перейти к вычислению искомых значений:

<math>\sigma_i = ||G(:,i)||_2</math>;

<math>U = [u_1,\dots,u_n]</math>, где <math>u_i = G(:,i)/\sigma_i</math>;

<math>V = J</math>, где <math>J</math> - накопленное произведение вращений Якоби.

Одностороннее вращение в плоскости <math>(j,k)</math> осуществляет процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''':

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

Здесь <math>R(j,k,\theta)</math> - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 244-245)</ref>.

                                                                                
<math>
\begin{matrix}
& & & & j & & k & & & \\
\end{matrix}
</math>

<math>
J_i = 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>a_{jk}</math> не слишком мал"''. Подходящим критерием останова является условие
<math>\begin{align}
|a_{jk}| \ge \varepsilon\sqrt{a_{jj}a_{kk}}
\end{align}</math>.

=== Вычислительное ядро алгоритма ===

Вычислительное ядро составляют множественные вычисления элементов матрицы <math>a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk}</math> и <math>a_{kk} = (G^TG)_{kk}</math>, а также вычисление промежуточной матрицы <math>GR(j, k, \theta)</math> путем умножения исходной матрицы на матрицу вращения Якоби в процессе выполнения процедуры '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''.

После приведения матрицы <math>G^TG</math> к виду, достаточно близкому к диагональному, вычисляются сингулярные числа:

<math>\sigma_i = ||G(:,i)||_2</math>

Матрица левых сингулярных векторов <math>U = [u_1,..., u_n]</math>:

<math>u_i = G(:,i)/\sigma_i</math>

И (если требуется) матрица правых сингулярных векторов <math>V</math> путем накопления произведения вращений Якоби:

<math>
V = J,
J = JR(j, k, \theta)
</math>

=== Макроструктура алгоритма ===
Основную часть метода составляет процедура односторонних вращений матрицы (процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''). Рассмотрим подробнее структуру этой процедуры.

- Вычисление <math>a_{jk}, a_{jj}, a_{kk} (a_{ij} = (G^TG)_{ij})</math>. В связи с тем, что нет необходимости каждый раз вычислять произведение <math>G^TG</math>, для вычисления каждого элемента требуется <math>n</math> операций умножения и <math>(n-1)</math> операций сложения.

- Вычисление <math>\tau, c, s, t</math> не представляют значительной вычислительной сложности.

- Одностороннее вращение матрицы (умножение входной матрицы на матрицу вращения) требует порядка <math>O(n)</math> операций.

В итоге, процедура '''''One-Sided-Jacobi-Rotation''''' требует <math>O(n)</math> операций.

После выполнения алгоритма односторонних вращений необходимо вычислить искомые матрицы <math>S, U, V</math> (матрица сингулярных чисел и матрицы левых и правых векторов).

- Вычисление сингулярных чисел <math>\sigma_i</math> требует <math>n</math> операций умножения.

- Вычисление <math>U</math> требует <math>O(n^2)</math> операций. Аналогично для матрицы <math>V</math>.

=== Схема реализации последовательного алгоритма ===
Метод можно описать следующим образом<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-263)</ref>:

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> - накопленное произведение вращений Якоби.

Метод Якоби для нахождения сингулярных чисел и векторов матрицы <math>G</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

=== Последовательная сложность алгоритма ===
Для вычисления сингулярных чисел и векторов матрицы порядка <math>n</math> в последовательном варианте на каждой итерации алгоритма вызывается процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''', сложность которой равна <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>j</math> — меняется в диапазоне от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения;
* <math>k</math> — меняется в диапазоне от <math>j+1</math> до <math>n</math>, принимая все целочисленные значения.

'''Вторая''' группа вершин соответствует вычислению значений <math>c</math> и <math>s</math> для матрицы вращения

'''Третья''' группа вершин соответствует вычислению произведений матриц <math>GR(k,j, \theta)</math> и <math>JR(k,j, \theta)</math>

<center>
<div class="thumb">
<div class="thumbinner" style="width:{{#expr: 3 * (400 + 35) + 4 * (3 - 1) + 8}}px">
<gallery widths=400px heights=400px>
File:JacobiSVG.png|Рисунок 1. Граф одной итерации алгоритма без отображения входных и выходных данных.
File:Akj.png|Рисунок 2. Граф вычисления значений элементов <math>a_{jj}, a_{jk}, a_{kk}</math>.
</gallery>
<gallery widths=400px heights=1200px>
File:JacobiSVD2.png|Рисунок 3. Граф одной итерации алгоритма с вычислением матрицы вращения <math>R</math> отображением входных и выходных данных.
</gallery>
<gallery widths=400px heights=400px>
File:Gij.png|Рисунок 4. Граф вычисления значений элементов матрицы <math>GR</math>
File:Jij.png|Рисунок 5. Граф вычисления значений элементов матрицы <math>JR</math>
</gallery>
</div>
</div>
</center>

=== Ресурс параллелизма алгоритма ===
В данном алгоритме существует возможность параллельного вычисления вращений Якоби<ref>Halil Snopce, Ilir Spahiu Parallel Computation of the SVD, pp. 458-460, http://www.wseas.us/e-library/conferences/2011/Paris/ECC/ECC-73.pdf</ref><ref>V. Strumpen, H. Hoffmann, A. Agarwal A Stream Algorithm for the SVD, pp. 6-7, 23-24, http://publications.csail.mit.edu/lcs/pubs/pdf/MIT-LCS-TM-641.pdf</ref>. Максимальное количество вращений в данном методе равно <math>\frac{n(n-1)}{2}</math>. В параллельной реализации, если <math>n</math> - четное, можно выполнить <math>\frac{n}{2}</math> вращений параллельно, если <math>n</math> - нечетное, то <math>\frac{n-1}{2}</math> вращений. Также следует учесть, что сложность самой процедуры вращения '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''' равна <math>O(n)</math>.

Таким образом, при классификации по высоте ЯПФ метод Якоби относится к алгоритмам со сложностью <math>O(n^2)</math>, а при классификации по ширине ЯПФ — к алгоритмам со сложностью <math>O(n)</math>.

=== Входные и выходные данные алгоритма ===
'''Входные данные''': плотная матрица <math>G</math> (элементы <math>g_{ij}</math>).
Дополнительные ограничения:
* <math>A= G^TG</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

'''Объём входных данных''': <math>n^2</math>

'''Выходные данные''': матрица <math>\Sigma = diag(\sigma_{i})</math>, где <math>\sigma_{i}</math> - сингулярные числа, матрица <math>U</math> левых сингулярных векторов и матрица <math>V</math> правых сингулярных векторов.

'''Объём выходных данных''': <math>2n^2 + n</math> (так как матрица сингулярных чисел <math>\Sigma</math> - диагональная, то достаточно хранить вектор сингулярных чисел (элементы <math>\sigma_i</math>), а также две матрицы <math>U, V</math> левых и правых сингулярных векторов).

=== Свойства алгоритма ===
Метод Якоби является самым медленным из имеющихся алгоритмов вычисления сингулярных чисел и сингулярных векторов матрицы. Тем не менее, для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. Кроме того, в методе Якоби заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, за счет этого метод достигает высокой относительной точности.

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов для метода Якоби является линейным, вычислительная мощность алгоритма также линейна.

Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности (число итераций зависит от входных данных и порогового значения).

== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===

=== Локальность данных и вычислений ===

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

=== Масштабируемость алгоритма и его реализации ===

==== Масштабируемость алгоритма ====
Оценки предела масштабируемости алгоритма были рассмотрены в разделе [[#Ресурс параллелизма алгоритма|Ресурс параллелизма алгоритма]].

==== Масштабируемость реализации алгоритма ====
Проведём исследование масштабируемости параллельной реализации метода Якоби нахождения сингулярных чисел и векторов. Исследование проводилось на суперкомпьютере "Ломоносов" [http://parallel.ru/cluster Суперкомпьютерного комплекса Московского университета].

Набор и границы значений изменяемых реализации алгоритма:

* число процессоров [1:16] с шагом 1;
* размер матрицы [100 : 900] с шагом 100.

На следующем рисунке приведен график среднего времени работы алгоритма в зависимости от изменяемых параметров запуска.

[[file: SQU3CTh5bI41.jpg | Рисунок 6. Параллельная реализация метода Якоби вычисления сингулярных чисел и векторов. Изменение времени выполнения в зависимости от числа процессоров и размера матрицы]]

При увеличении числа процессов эффективность на рассмотренной области изменений параметров запуска падает.

[https://software.intel.com/ru-ru/node/521152#248AD4D8-005A-4800-9A73-428952BC1E44 Исследованная параллельная реализация, реализованная в библиотеке Intel MKL]

=== Динамические характеристики и эффективность реализации алгоритма ===

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

=== Существующие реализации алгоритма ===
Несмотря на то, что метод Якоби является самым медленным из методов вычисления сингулярных чисел и векторов матрицы, он реализован в некоторых пакетах, например, LAPACK, ScaLAPACK, Intel MKL<ref>Intel® Math Kernel Library, https://software.intel.com/ru-ru/node/521153#CAF8CA08-1C69-48F8-B0C3-B2B27AFA92FA</ref>.

== Литература ==

Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов

2016-10-29T09:28:16Z

Светлана:

{{Assignment}}
{{algorithm
| name = Метод Якоби вычисления сингулярных чисел и векторов
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math>O(n^2)</math>
| pf_width = <math>O(n)</math>
| input_data = <math>n^2</math>
| output_data = <math>2n^2 + n</math>
}}
Авторы статьи: [[Участник:Почернина Елена|Почернина Елена]], группа 601 (разделы [[#Общее описание алгоритма|1.1]], [[#Математическое описание алгоритма|1.2]], [[#Макроструктура алгоритма|1.4]], [[#Схема реализации последовательного алгоритма|1.5]], [[#Последовательная сложность алгоритма|1.6]], [[#Ресурс параллелизма алгоритма|1.8]], [[#Существующие реализации алгоритма|2.7]]), [[Участница:Светлана|Костюкова Светлана]], группа 601 (разделы [[#Вычислительное ядро алгоритма|1.3]], [[#Схема реализации последовательного алгоритма|1.5]], [[#Информационный граф|1.7]], [[#Входные и выходные данные алгоритма|1.9]], [[#Свойства алгоритма|1.10]])

== Свойства и структура алгоритма ==

=== Общее описание алгоритма ===

'''Метод Якоби''' позволяет найти сингулярное разложение плотной матрицы <math>G</math> размера <math>n</math>×<math>n</math>. Такое разложение можно записать в виде <math>G = U\Sigma V^*</math>, где <math>U, V</math> - унитарные (''в вещественном случае - ортогональные'') матрицы размера <math>n</math>×<math>n</math>, <math>\Sigma</math> - диагональная матрица размера <math>n</math>×<math>n</math> с вещественными положительными числами на главной диагонали. Диагональные элементы матрицы <math>\Sigma</math> называются сингулярными числами матрицы <math>G</math>, а столбцы матриц <math>U,V</math> - левыми и правыми сингулярными векторами.

''Неявный''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref> метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице <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>. Поэтому алгоритм называется ''методом односторонних вращений''.

Метод Якоби для вычисления сингулярных чисел и векторов является самым медленным из имеющихся алгоритмов поиска сингулярных чисел и векторов, но тем не менее, интерес к нему сохраняется. Для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. В частности, метод Якоби вычисляет сингулярные числа матрицы <math>G</math> с высокой точностью, если G может быть представлена в виде <math>G = DX</math>, где <math>D</math> – диагональная матрица, а <math>X</math> – хорошо обусловлена. В этом случае заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, в то время как другие алгоритмы включают в себя приведение матрицы к двухдиагональной форме, из-за чего и теряют все верные разряды во всех сингулярных числах, кроме старшего.

=== Математическое описание алгоритма ===
Исходные данные: матрица <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> - матрица правых сингулярных векторов.

Основной алгоритм заключается в следующем<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref>:

Для всех элементов матрицы <math>A^{(i)} = (G^TG)^{(i)}</math>, находящихся вне главной диагонали <math>(a_{12}, a_{13},\dots, a_{1n}, a_{23}, a_{24},\dots, a_{2n},\dots, a_{n-1n})</math> проверяем условие <math>|A_{jk}|<\varepsilon</math>:
Если <math>|A_{jk}^{(i)}| \le \varepsilon</math>, то переходим к следующему элементу,
Если <math>|A_{jk}^{(i)}| > \varepsilon</math>, то нужно выполнить вращение в плоскости <math>(j,k)</math>, после чего перейти к следующей итерации алгоритма.
Если на очередной итерации не было совершено ни одного вращения, значит матрица <math>A</math> стала достаточна близка к диагональной, и можно перейти к вычислению искомых значений:

<math>\sigma_i = ||G(:,i)||_2</math>;

<math>U = [u_1,\dots,u_n]</math>, где <math>u_i = G(:,i)/\sigma_i</math>;

<math>V = J</math>, где <math>J</math> - накопленное произведение вращений Якоби.

Одностороннее вращение в плоскости <math>(j,k)</math> осуществляет процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''':

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

Здесь <math>R(j,k,\theta)</math> - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 244-245)</ref>.

                                                                                
<math>
\begin{matrix}
& & & & j & & k & & & \\
\end{matrix}
</math>

<math>
J_i = 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>a_{jk}</math> не слишком мал"''. Подходящим критерием останова является условие
<math>\begin{align}
|a_{jk}| \ge \varepsilon\sqrt{a_{jj}a_{kk}}
\end{align}</math>.

=== Вычислительное ядро алгоритма ===

Вычислительное ядро составляют множественные вычисления элементов матрицы <math>a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk}</math> и <math>a_{kk} = (G^TG)_{kk}</math>, а также вычисление промежуточной матрицы <math>GR(j, k, \theta)</math> путем умножения исходной матрицы на матрицу вращения Якоби в процессе выполнения процедуры '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''.

После приведения матрицы <math>G^TG</math> к виду, достаточно близкому к диагональному, вычисляются сингулярные числа:

<math>\sigma_i = ||G(:,i)||_2</math>

Матрица левых сингулярных векторов <math>U = [u_1,..., u_n]</math>:

<math>u_i = G(:,i)/\sigma_i</math>

И (если требуется) матрица правых сингулярных векторов <math>V</math> путем накопления произведения вращений Якоби:

<math>
V = J,
J = JR(j, k, \theta)
</math>

=== Макроструктура алгоритма ===
Основную часть метода составляет процедура односторонних вращений матрицы (процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''). Рассмотрим подробнее структуру этой процедуры.

- Вычисление <math>a_{jk}, a_{jj}, a_{kk} (a_{ij} = (G^TG)_{ij})</math>. В связи с тем, что нет необходимости каждый раз вычислять произведение <math>G^TG</math>, для вычисления каждого элемента требуется <math>n</math> операций умножения и <math>(n-1)</math> операций сложения.

- Вычисление <math>\tau, c, s, t</math> не представляют значительной вычислительной сложности.

- Одностороннее вращение матрицы (умножение входной матрицы на матрицу вращения) требует порядка <math>O(n)</math> операций.

В итоге, процедура '''''One-Sided-Jacobi-Rotation''''' требует <math>O(n)</math> операций.

После выполнения алгоритма односторонних вращений необходимо вычислить искомые матрицы <math>S, U, V</math> (матрица сингулярных чисел и матрицы левых и правых векторов).

- Вычисление сингулярных чисел <math>\sigma_i</math> требует <math>n</math> операций умножения.

- Вычисление <math>U</math> требует <math>O(n^2)</math> операций. Аналогично для матрицы <math>V</math>.

=== Схема реализации последовательного алгоритма ===
Метод можно описать следующим образом<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-263)</ref>:

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> - накопленное произведение вращений Якоби.

Метод Якоби для нахождения сингулярных чисел и векторов матрицы <math>G</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

=== Последовательная сложность алгоритма ===
Для вычисления сингулярных чисел и векторов матрицы порядка <math>n</math> в последовательном варианте на каждой итерации алгоритма вызывается процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''', сложность которой равна <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>j</math> — меняется в диапазоне от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения;
* <math>k</math> — меняется в диапазоне от <math>j+1</math> до <math>n</math>, принимая все целочисленные значения.

'''Вторая''' группа вершин соответствует вычислению значений <math>c</math> и <math>s</math> для матрицы вращения

'''Третья''' группа вершин соответствует вычислению произведений матриц <math>GR(k,j, \theta)</math> и <math>JR(k,j, \theta)</math>

<center>
<div class="thumb">
<div class="thumbinner" style="width:{{#expr: 3 * (400 + 35) + 4 * (3 - 1) + 8}}px">
<gallery widths=400px heights=400px>
File:JacobiSVG.png|Рисунок 1. Граф одной итерации алгоритма без отображения входных и выходных данных.
File:Akj.png|Рисунок 2. Граф вычисления значений элементов <math>a_{jj}, a_{jk}, a_{kk}</math>.
</gallery>
<gallery widths=400px heights=1200px>
File:JacobiSVD2.png|Рисунок 3. Граф одной итерации алгоритма с вычислением матрицы вращения <math>R</math> отображением входных и выходных данных.
</gallery>
<gallery widths=400px heights=400px>
File:Gij.png|Рисунок 4. Граф вычисления значений элементов матрицы <math>GR</math>
File:Jij.png|Рисунок 5. Граф вычисления значений элементов матрицы <math>JR</math>
</galery>
</div>
</div>
</center>

=== Ресурс параллелизма алгоритма ===
В данном алгоритме существует возможность параллельного вычисления вращений Якоби<ref>Halil Snopce, Ilir Spahiu Parallel Computation of the SVD, pp. 458-460, http://www.wseas.us/e-library/conferences/2011/Paris/ECC/ECC-73.pdf</ref><ref>V. Strumpen, H. Hoffmann, A. Agarwal A Stream Algorithm for the SVD, pp. 6-7, 23-24, http://publications.csail.mit.edu/lcs/pubs/pdf/MIT-LCS-TM-641.pdf</ref>. Максимальное количество вращений в данном методе равно <math>\frac{n(n-1)}{2}</math>. В параллельной реализации, если <math>n</math> - четное, можно выполнить <math>\frac{n}{2}</math> вращений параллельно, если <math>n</math> - нечетное, то <math>\frac{n-1}{2}</math> вращений. Также следует учесть, что сложность самой процедуры вращения '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''' равна <math>O(n)</math>.

Таким образом, при классификации по высоте ЯПФ метод Якоби относится к алгоритмам со сложностью <math>O(n^2)</math>, а при классификации по ширине ЯПФ — к алгоритмам со сложностью <math>O(n)</math>.

=== Входные и выходные данные алгоритма ===
'''Входные данные''': плотная матрица <math>G</math> (элементы <math>g_{ij}</math>).
Дополнительные ограничения:
* <math>A= G^TG</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

'''Объём входных данных''': <math>n^2</math>

'''Выходные данные''': матрица <math>\Sigma = diag(\sigma_{i})</math>, где <math>\sigma_{i}</math> - сингулярные числа, матрица <math>U</math> левых сингулярных векторов и матрица <math>V</math> правых сингулярных векторов.

'''Объём выходных данных''': <math>2n^2 + n</math> (так как матрица сингулярных чисел <math>\Sigma</math> - диагональная, то достаточно хранить вектор сингулярных чисел (элементы <math>\sigma_i</math>), а также две матрицы <math>U, V</math> левых и правых сингулярных векторов).

=== Свойства алгоритма ===
Метод Якоби является самым медленным из имеющихся алгоритмов вычисления сингулярных чисел и сингулярных векторов матрицы. Тем не менее, для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. Кроме того, в методе Якоби заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, за счет этого метод достигает высокой относительной точности.

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов для метода Якоби является линейным, вычислительная мощность алгоритма также линейна.

Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности (число итераций зависит от входных данных и порогового значения).

== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===

=== Локальность данных и вычислений ===

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

=== Масштабируемость алгоритма и его реализации ===
Данный раздел будет дополнен позже в связи с продлением срока сдачи (до 15 ноября).

=== Динамические характеристики и эффективность реализации алгоритма ===

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

=== Существующие реализации алгоритма ===
Несмотря на то, что метод Якоби является самым медленным из методов вычисления сингулярных чисел и векторов матрицы, он реализован в некоторых пакетах, например, LAPACK, ScaLAPACK, Intel MKL<ref>Intel® Math Kernel Library, https://software.intel.com/ru-ru/node/521153#CAF8CA08-1C69-48F8-B0C3-B2B27AFA92FA</ref>.

== Литература ==

Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов

2016-10-29T09:22:52Z

Светлана: /* Информационный граф */

{{algorithm
| name = Метод Якоби вычисления сингулярных чисел и векторов
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math>O(n^2)</math>
| pf_width = <math>O(n)</math>
| input_data = <math>n^2</math>
| output_data = <math>2n^2 + n</math>
}}
Авторы статьи: [[Участник:Почернина Елена|Почернина Елена]], группа 601 (разделы [[#Общее описание алгоритма|1.1]], [[#Математическое описание алгоритма|1.2]], [[#Макроструктура алгоритма|1.4]], [[#Схема реализации последовательного алгоритма|1.5]], [[#Последовательная сложность алгоритма|1.6]], [[#Ресурс параллелизма алгоритма|1.8]], [[#Существующие реализации алгоритма|2.7]]), [[Участница:Светлана|Костюкова Светлана]], группа 601 (разделы [[#Вычислительное ядро алгоритма|1.3]], [[#Схема реализации последовательного алгоритма|1.5]], [[#Информационный граф|1.7]], [[#Входные и выходные данные алгоритма|1.9]], [[#Свойства алгоритма|1.10]])

== Свойства и структура алгоритма ==

=== Общее описание алгоритма ===

'''Метод Якоби''' позволяет найти сингулярное разложение плотной матрицы <math>G</math> размера <math>n</math>×<math>n</math>. Такое разложение можно записать в виде <math>G = U\Sigma V^*</math>, где <math>U, V</math> - унитарные (''в вещественном случае - ортогональные'') матрицы размера <math>n</math>×<math>n</math>, <math>\Sigma</math> - диагональная матрица размера <math>n</math>×<math>n</math> с вещественными положительными числами на главной диагонали. Диагональные элементы матрицы <math>\Sigma</math> называются сингулярными числами матрицы <math>G</math>, а столбцы матриц <math>U,V</math> - левыми и правыми сингулярными векторами.

''Неявный''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref> метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице <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>. Поэтому алгоритм называется ''методом односторонних вращений''.

Метод Якоби для вычисления сингулярных чисел и векторов является самым медленным из имеющихся алгоритмов поиска сингулярных чисел и векторов, но тем не менее, интерес к нему сохраняется. Для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. В частности, метод Якоби вычисляет сингулярные числа матрицы <math>G</math> с высокой точностью, если G может быть представлена в виде <math>G = DX</math>, где <math>D</math> – диагональная матрица, а <math>X</math> – хорошо обусловлена. В этом случае заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, в то время как другие алгоритмы включают в себя приведение матрицы к двухдиагональной форме, из-за чего и теряют все верные разряды во всех сингулярных числах, кроме старшего.

=== Математическое описание алгоритма ===
Исходные данные: матрица <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> - матрица правых сингулярных векторов.

Основной алгоритм заключается в следующем<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref>:

Для всех элементов матрицы <math>A^{(i)} = (G^TG)^{(i)}</math>, находящихся вне главной диагонали <math>(a_{12}, a_{13},\dots, a_{1n}, a_{23}, a_{24},\dots, a_{2n},\dots, a_{n-1n})</math> проверяем условие <math>|A_{jk}|<\varepsilon</math>:
Если <math>|A_{jk}^{(i)}| \le \varepsilon</math>, то переходим к следующему элементу,
Если <math>|A_{jk}^{(i)}| > \varepsilon</math>, то нужно выполнить вращение в плоскости <math>(j,k)</math>, после чего перейти к следующей итерации алгоритма.
Если на очередной итерации не было совершено ни одного вращения, значит матрица <math>A</math> стала достаточна близка к диагональной, и можно перейти к вычислению искомых значений:

<math>\sigma_i = ||G(:,i)||_2</math>;

<math>U = [u_1,\dots,u_n]</math>, где <math>u_i = G(:,i)/\sigma_i</math>;

<math>V = J</math>, где <math>J</math> - накопленное произведение вращений Якоби.

Одностороннее вращение в плоскости <math>(j,k)</math> осуществляет процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''':

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

Здесь <math>R(j,k,\theta)</math> - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 244-245)</ref>.

                                                                                
<math>
\begin{matrix}
& & & & j & & k & & & \\
\end{matrix}
</math>

<math>
J_i = 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>a_{jk}</math> не слишком мал"''. Подходящим критерием останова является условие
<math>\begin{align}
|a_{jk}| \ge \varepsilon\sqrt{a_{jj}a_{kk}}
\end{align}</math>.

=== Вычислительное ядро алгоритма ===

Вычислительное ядро составляют множественные вычисления элементов матрицы <math>a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk}</math> и <math>a_{kk} = (G^TG)_{kk}</math>, а также вычисление промежуточной матрицы <math>GR(j, k, \theta)</math> путем умножения исходной матрицы на матрицу вращения Якоби в процессе выполнения процедуры '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''.

После приведения матрицы <math>G^TG</math> к виду, достаточно близкому к диагональному, вычисляются сингулярные числа:

<math>\sigma_i = ||G(:,i)||_2</math>

Матрица левых сингулярных векторов <math>U = [u_1,..., u_n]</math>:

<math>u_i = G(:,i)/\sigma_i</math>

И (если требуется) матрица правых сингулярных векторов <math>V</math> путем накопления произведения вращений Якоби:

<math>
V = J,
J = JR(j, k, \theta)
</math>

=== Макроструктура алгоритма ===
Основную часть метода составляет процедура односторонних вращений матрицы (процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''). Рассмотрим подробнее структуру этой процедуры.

- Вычисление <math>a_{jk}, a_{jj}, a_{kk} (a_{ij} = (G^TG)_{ij})</math>. В связи с тем, что нет необходимости каждый раз вычислять произведение <math>G^TG</math>, для вычисления каждого элемента требуется <math>n</math> операций умножения и <math>(n-1)</math> операций сложения.

- Вычисление <math>\tau, c, s, t</math> не представляют значительной вычислительной сложности.

- Одностороннее вращение матрицы (умножение входной матрицы на матрицу вращения) требует порядка <math>O(n)</math> операций.

В итоге, процедура '''''One-Sided-Jacobi-Rotation''''' требует <math>O(n)</math> операций.

После выполнения алгоритма односторонних вращений необходимо вычислить искомые матрицы <math>S, U, V</math> (матрица сингулярных чисел и матрицы левых и правых векторов).

- Вычисление сингулярных чисел <math>\sigma_i</math> требует <math>n</math> операций умножения.

- Вычисление <math>U</math> требует <math>O(n^2)</math> операций. Аналогично для матрицы <math>V</math>.

=== Схема реализации последовательного алгоритма ===
Метод можно описать следующим образом<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-263)</ref>:

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> - накопленное произведение вращений Якоби.

Метод Якоби для нахождения сингулярных чисел и векторов матрицы <math>G</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

=== Последовательная сложность алгоритма ===
Для вычисления сингулярных чисел и векторов матрицы порядка <math>n</math> в последовательном варианте на каждой итерации алгоритма вызывается процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''', сложность которой равна <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>j</math> — меняется в диапазоне от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения;
* <math>k</math> — меняется в диапазоне от <math>j+1</math> до <math>n</math>, принимая все целочисленные значения.

'''Вторая''' группа вершин соответствует вычислению значений <math>c</math> и <math>s</math> для матрицы вращения

'''Третья''' группа вершин соответствует вычислению произведений матриц <math>GR(k,j, \theta)</math> и <math>JR(k,j, \theta)</math>

<center>
<div class="thumb">
<div class="thumbinner" style="width:{{#expr: 3 * (400 + 35) + 4 * (3 - 1) + 8}}px">
<gallery widths=400px heights=400px>
File:JacobiSVG.png|Рисунок 1. Граф одной итерации алгоритма без отображения входных и выходных данных.
File:Akj.png|Рисунок 2. Граф вычисления значений элементов <math>a_{jj}, a_{jk}, a_{kk}</math>.
</gallery>
<gallery widths=400px heights=1200px>
File:JacobiSVD2.png|Рисунок 3. Граф одной итерации алгоритма с вычислением матрицы вращения <math>R</math> отображением входных и выходных данных.
</gallery>
<gallery widths=400px heights=400px>
File:Gij.png|Рисунок 4. Граф вычисления значений элементов матрицы <math>GR</math>
File:Jij.png|Рисунок 5. Граф вычисления значений элементов матрицы <math>JR</math>
</galery>
</div>
</div>
</center>

=== Ресурс параллелизма алгоритма ===
В данном алгоритме существует возможность параллельного вычисления вращений Якоби<ref>Halil Snopce, Ilir Spahiu Parallel Computation of the SVD, pp. 458-460, http://www.wseas.us/e-library/conferences/2011/Paris/ECC/ECC-73.pdf</ref><ref>V. Strumpen, H. Hoffmann, A. Agarwal A Stream Algorithm for the SVD, pp. 6-7, 23-24, http://publications.csail.mit.edu/lcs/pubs/pdf/MIT-LCS-TM-641.pdf</ref>. Максимальное количество вращений в данном методе равно <math>\frac{n(n-1)}{2}</math>. В параллельной реализации, если <math>n</math> - четное, можно выполнить <math>\frac{n}{2}</math> вращений параллельно, если <math>n</math> - нечетное, то <math>\frac{n-1}{2}</math> вращений. Также следует учесть, что сложность самой процедуры вращения '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''' равна <math>O(n)</math>.

Таким образом, при классификации по высоте ЯПФ метод Якоби относится к алгоритмам со сложностью <math>O(n^2)</math>, а при классификации по ширине ЯПФ — к алгоритмам со сложностью <math>O(n)</math>.

=== Входные и выходные данные алгоритма ===
'''Входные данные''': плотная матрица <math>G</math> (элементы <math>g_{ij}</math>).
Дополнительные ограничения:
* <math>A= G^TG</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

'''Объём входных данных''': <math>n^2</math>

'''Выходные данные''': матрица <math>\Sigma = diag(\sigma_{i})</math>, где <math>\sigma_{i}</math> - сингулярные числа, матрица <math>U</math> левых сингулярных векторов и матрица <math>V</math> правых сингулярных векторов.

'''Объём выходных данных''': <math>2n^2 + n</math> (так как матрица сингулярных чисел <math>\Sigma</math> - диагональная, то достаточно хранить вектор сингулярных чисел (элементы <math>\sigma_i</math>), а также две матрицы <math>U, V</math> левых и правых сингулярных векторов).

=== Свойства алгоритма ===
Метод Якоби является самым медленным из имеющихся алгоритмов вычисления сингулярных чисел и сингулярных векторов матрицы. Тем не менее, для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. Кроме того, в методе Якоби заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, за счет этого метод достигает высокой относительной точности.

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов для метода Якоби является линейным, вычислительная мощность алгоритма также линейна.

Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности (число итераций зависит от входных данных и порогового значения).

== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===

=== Локальность данных и вычислений ===

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

=== Масштабируемость алгоритма и его реализации ===
Данный раздел будет дополнен позже в связи с продлением срока сдачи (до 15 ноября).

=== Динамические характеристики и эффективность реализации алгоритма ===

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

=== Существующие реализации алгоритма ===
Несмотря на то, что метод Якоби является самым медленным из методов вычисления сингулярных чисел и векторов матрицы, он реализован в некоторых пакетах, например, LAPACK, ScaLAPACK, Intel MKL<ref>Intel® Math Kernel Library, https://software.intel.com/ru-ru/node/521153#CAF8CA08-1C69-48F8-B0C3-B2B27AFA92FA</ref>.

== Литература ==

Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов

2016-10-27T15:41:44Z

Светлана: /* Информационный граф */

{{algorithm
| name = Метод Якоби вычисления сингулярных чисел и векторов
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math>O(n^2)</math>
| pf_width = <math>O(n)</math>
| input_data = <math>n^2</math>
| output_data = <math>2n^2 + n</math>
}}
Авторы статьи: [[Участник:Почернина Елена|Почернина Елена]], группа 601 (разделы [[#Общее описание алгоритма|1.1]], [[#Математическое описание алгоритма|1.2]], [[#Макроструктура алгоритма|1.4]], [[#Схема реализации последовательного алгоритма|1.5]], [[#Последовательная сложность алгоритма|1.6]], [[#Ресурс параллелизма алгоритма|1.8]], [[#Существующие реализации алгоритма|2.7]]), [[Участница:Светлана|Костюкова Светлана]], группа 601 (разделы [[#Вычислительное ядро алгоритма|1.3]], [[#Схема реализации последовательного алгоритма|1.5]], [[#Информационный граф|1.7]], [[#Входные и выходные данные алгоритма|1.9]], [[#Свойства алгоритма|1.10]])

== Свойства и структура алгоритма ==

=== Общее описание алгоритма ===

'''Метод Якоби''' позволяет найти сингулярное разложение плотной матрицы <math>G</math> размера <math>n</math>×<math>n</math>. Такое разложение можно записать в виде <math>G = U\Sigma V^*</math>, где <math>U, V</math> - унитарные (''в вещественном случае - ортогональные'') матрицы размера <math>n</math>×<math>n</math>, <math>\Sigma</math> - диагональная матрица размера <math>n</math>×<math>n</math> с вещественными положительными числами на главной диагонали. Диагональные элементы матрицы <math>\Sigma</math> называются сингулярными числами матрицы <math>G</math>, а столбцы матриц <math>U,V</math> - левыми и правыми сингулярными векторами.

''Неявный''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref> метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице <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>. Поэтому алгоритм называется ''методом односторонних вращений''.

Метод Якоби для вычисления сингулярных чисел и векторов является самым медленным из имеющихся алгоритмов поиска сингулярных чисел и векторов, но тем не менее, интерес к нему сохраняется. Для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. В частности, метод Якоби вычисляет сингулярные числа матрицы <math>G</math> с высокой точностью, если G может быть представлена в виде <math>G = DX</math>, где <math>D</math> – диагональная матрица, а <math>X</math> – хорошо обусловлена. В этом случае заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, в то время как другие алгоритмы включают в себя приведение матрицы к двухдиагональной форме, из-за чего и теряют все верные разряды во всех сингулярных числах, кроме старшего.

=== Математическое описание алгоритма ===
Исходные данные: матрица <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> - матрица правых сингулярных векторов.

Основной алгоритм заключается в следующем<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref>:

Для всех элементов матрицы <math>A^{(i)} = (G^TG)^{(i)}</math>, находящихся вне главной диагонали <math>(a_{12}, a_{13},\dots, a_{1n}, a_{23}, a_{24},\dots, a_{2n},\dots, a_{n-1n})</math> проверяем условие <math>|A_{jk}|<\varepsilon</math>:
Если <math>|A_{jk}^{(i)}| \le \varepsilon</math>, то переходим к следующему элементу,
Если <math>|A_{jk}^{(i)}| > \varepsilon</math>, то нужно выполнить вращение в плоскости <math>(j,k)</math>, после чего перейти к следующей итерации алгоритма.
Если на очередной итерации не было совершено ни одного вращения, значит матрица <math>A</math> стала достаточна близка к диагональной, и можно перейти к вычислению искомых значений:

<math>\sigma_i = ||G(:,i)||_2</math>;

<math>U = [u_1,\dots,u_n]</math>, где <math>u_i = G(:,i)/\sigma_i</math>;

<math>V = J</math>, где <math>J</math> - накопленное произведение вращений Якоби.

Одностороннее вращение в плоскости <math>(j,k)</math> осуществляет процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''':

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

Здесь <math>R(j,k,\theta)</math> - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 244-245)</ref>.

                                                                                
<math>
\begin{matrix}
& & & & j & & k & & & \\
\end{matrix}
</math>

<math>
J_i = 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>a_{jk}</math> не слишком мал"''. Подходящим критерием останова является условие
<math>\begin{align}
|a_{jk}| \ge \varepsilon\sqrt{a_{jj}a_{kk}}
\end{align}</math>.

=== Вычислительное ядро алгоритма ===

Вычислительное ядро составляют множественные вычисления элементов матрицы <math>a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk}</math> и <math>a_{kk} = (G^TG)_{kk}</math>, а также вычисление промежуточной матрицы <math>GR(j, k, \theta)</math> путем умножения исходной матрицы на матрицу вращения Якоби в процессе выполнения процедуры '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''.

После приведения матрицы <math>G^TG</math> к виду, достаточно близкому к диагональному, вычисляются сингулярные числа:

<math>\sigma_i = ||G(:,i)||_2</math>

Матрица левых сингулярных векторов <math>U = [u_1,..., u_n]</math>:

<math>u_i = G(:,i)/\sigma_i</math>

И (если требуется) матрица правых сингулярных векторов <math>V</math> путем накопления произведения вращений Якоби:

<math>
V = J,
J = JR(j, k, \theta)
</math>

=== Макроструктура алгоритма ===
Основную часть метода составляет процедура односторонних вращений матрицы (процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''). Рассмотрим подробнее структуру этой процедуры.

- Вычисление <math>a_{jk}, a_{jj}, a_{kk} (a_{ij} = (G^TG)_{ij})</math>. В связи с тем, что нет необходимости каждый раз вычислять произведение <math>G^TG</math>, для вычисления каждого элемента требуется <math>n</math> операций умножения и <math>(n-1)</math> операций сложения.

- Вычисление <math>\tau, c, s, t</math> не представляют значительной вычислительной сложности.

- Одностороннее вращение матрицы (умножение входной матрицы на матрицу вращения) требует порядка <math>O(n)</math> операций.

В итоге, процедура '''''One-Sided-Jacobi-Rotation''''' требует <math>O(n)</math> операций.

После выполнения алгоритма односторонних вращений необходимо вычислить искомые матрицы <math>S, U, V</math> (матрица сингулярных чисел и матрицы левых и правых векторов).

- Вычисление сингулярных чисел <math>\sigma_i</math> требует <math>n</math> операций умножения.

- Вычисление <math>U</math> требует <math>O(n^2)</math> операций. Аналогично для матрицы <math>V</math>.

=== Схема реализации последовательного алгоритма ===
Метод можно описать следующим образом<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-263)</ref>:

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> - накопленное произведение вращений Якоби.

Метод Якоби для нахождения сингулярных чисел и векторов матрицы <math>G</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

=== Последовательная сложность алгоритма ===
Для вычисления сингулярных чисел и векторов матрицы порядка <math>n</math> в последовательном варианте на каждой итерации алгоритма вызывается процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''', сложность которой равна <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>j</math> — меняется в диапазоне от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения;
* <math>k</math> — меняется в диапазоне от <math>j+1</math> до <math>n</math>, принимая все целочисленные значения.

'''Вторая''' группа вершин соответствует вычислению значений <math>c</math> и <math>s</math> для матрицы вращения

'''Третья''' группа вершин соответствует вычислению произведений матриц <math>GR(k,j, \theta)</math> и <math>JR(k,j, \theta)</math>

<center>
<div class="thumb">
<div class="thumbinner" style="width:{{#expr: 3 * (400 + 35) + 4 * (3 - 1) + 8}}px">
<gallery widths=400px heights=400px>
File:JacobiSVG.png|Рисунок 1. Граф одной итерации алгоритма без отображения входных и выходных данных.
File:Akj.png|Рисунок 2. Граф вычисления значений элементов <math>a_{jj}, a_{jk}, a_{kk}</math>.
</gallery>
<gallery widths=400px heights=1200px>
File:JacobiSVD2.png|Рисунок 3. Граф одной итерации алгоритма с вычислением матрицы вращения <math>R</math> отображением входных и выходных данных.
</gallery>
</div>
</div>
</center>

=== Ресурс параллелизма алгоритма ===
В данном алгоритме существует возможность параллельного вычисления вращений Якоби<ref>Halil Snopce, Ilir Spahiu Parallel Computation of the SVD, pp. 458-460, http://www.wseas.us/e-library/conferences/2011/Paris/ECC/ECC-73.pdf</ref><ref>V. Strumpen, H. Hoffmann, A. Agarwal A Stream Algorithm for the SVD, pp. 6-7, 23-24, http://publications.csail.mit.edu/lcs/pubs/pdf/MIT-LCS-TM-641.pdf</ref>. Максимальное количество вращений в данном методе равно <math>\frac{n(n-1)}{2}</math>. В параллельной реализации, если <math>n</math> - четное, можно выполнить <math>\frac{n}{2}</math> вращений параллельно, если <math>n</math> - нечетное, то <math>\frac{n-1}{2}</math> вращений. Также следует учесть, что сложность самой процедуры вращения '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''' равна <math>O(n)</math>.

Таким образом, при классификации по высоте ЯПФ метод Якоби относится к алгоритмам со сложностью <math>O(n^2)</math>, а при классификации по ширине ЯПФ — к алгоритмам со сложностью <math>O(n)</math>.

=== Входные и выходные данные алгоритма ===
'''Входные данные''': плотная матрица <math>G</math> (элементы <math>g_{ij}</math>).
Дополнительные ограничения:
* <math>A= G^TG</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

'''Объём входных данных''': <math>n^2</math>

'''Выходные данные''': матрица <math>\Sigma = diag(\sigma_{i})</math>, где <math>\sigma_{i}</math> - сингулярные числа, матрица <math>U</math> левых сингулярных векторов и матрица <math>V</math> правых сингулярных векторов.

'''Объём выходных данных''': <math>2n^2 + n</math> (так как матрица сингулярных чисел <math>\Sigma</math> - диагональная, то достаточно хранить вектор сингулярных чисел (элементы <math>\sigma_i</math>), а также две матрицы <math>U, V</math> левых и правых сингулярных векторов).

=== Свойства алгоритма ===
Метод Якоби является самым медленным из имеющихся алгоритмов вычисления сингулярных чисел и сингулярных векторов матрицы. Тем не менее, для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. Кроме того, в методе Якоби заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, за счет этого метод достигает высокой относительной точности.

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов для метода Якоби является линейным, вычислительная мощность алгоритма также линейна.

Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности (число итераций зависит от входных данных и порогового значения).

== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===

=== Локальность данных и вычислений ===

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

=== Масштабируемость алгоритма и его реализации ===
Данный раздел будет дополнен позже в связи с продлением срока сдачи (до 15 ноября).

=== Динамические характеристики и эффективность реализации алгоритма ===

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

=== Существующие реализации алгоритма ===
Несмотря на то, что метод Якоби является самым медленным из методов вычисления сингулярных чисел и векторов матрицы, он реализован в некоторых пакетах, например, LAPACK, ScaLAPACK, Intel MKL<ref>Intel® Math Kernel Library, https://software.intel.com/ru-ru/node/521153#CAF8CA08-1C69-48F8-B0C3-B2B27AFA92FA</ref>.

== Литература ==

Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов

2016-10-27T14:35:16Z

Светлана: /* Информационный граф */

{{algorithm
| name = Метод Якоби вычисления сингулярных чисел и векторов
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math>O(n^2)</math>
| pf_width = <math>O(n)</math>
| input_data = <math>n^2</math>
| output_data = <math>2n^2 + n</math>
}}
Авторы статьи: [[Участник:Почернина Елена|Почернина Елена]] (группа 601), [[Участница:Светлана|Костюкова Светлана]] (группа 601)

== Свойства и структура алгоритма ==

=== Общее описание алгоритма ===

'''Метод Якоби''' позволяет найти сингулярное разложение плотной матрицы <math>G</math> размера <math>n</math>×<math>n</math>. Такое разложение можно записать в виде <math>G = U\Sigma V^*</math>, где <math>U, V</math> - унитарные (''в вещественном случае - ортогональные'') матрицы размера <math>n</math>×<math>n</math>, <math>\Sigma</math> - диагональная матрица размера <math>n</math>×<math>n</math> с вещественными положительными числами на главной диагонали. Диагональные элементы матрицы <math>\Sigma</math> называются сингулярными числами матрицы <math>G</math>, а столбцы матриц <math>U,V</math> - левыми и правыми сингулярными векторами.

''Неявный''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref> метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице <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>. Поэтому алгоритм называется ''методом односторонних вращений''.

Метод Якоби для вычисления сингулярных чисел и векторов является самым медленным из имеющихся алгоритмов поиска сингулярных чисел и векторов, но тем не менее, интерес к нему сохраняется. Для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. В частности, метод Якоби вычисляет сингулярные числа матрицы <math>G</math> с высокой точностью, если G может быть представлена в виде <math>G = DX</math>, где <math>D</math> – диагональная матрица, а <math>X</math> – хорошо обусловлена. В этом случае заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, в то время как другие алгоритмы включают в себя приведение матрицы к двухдиагональной форме, из-за чего и теряют все верные разряды во всех сингулярных числах, кроме старшего.

=== Математическое описание алгоритма ===
Исходные данные: матрица <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> - матрица правых сингулярных векторов.

Основной алгоритм заключается в следующем<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref>:

Для всех элементов матрицы <math>A^{(i)} = (G^TG)^{(i)}</math>, находящихся вне главной диагонали <math>(a_{12}, a_{13},\dots, a_{1n}, a_{23}, a_{24},\dots, a_{2n},\dots, a_{n-1n})</math> проверяем условие <math>|A_{jk}|<\varepsilon</math>:
Если <math>|A_{jk}^{(i)}| \le \varepsilon</math>, то переходим к следующему элементу,
Если <math>|A_{jk}^{(i)}| > \varepsilon</math>, то нужно выполнить вращение в плоскости <math>(j,k)</math>, после чего перейти к следующей итерации алгоритма.
Если на очередной итерации не было совершено ни одного вращения, значит матрица <math>A</math> стала достаточна близка к диагональной, и можно перейти к вычислению искомых значений:

<math>\sigma_i = ||G(:,i)||_2</math>;

<math>U = [u_1,\dots,u_n]</math>, где <math>u_i = G(:,i)/\sigma_i</math>;

<math>V = J</math>, где <math>J</math> - накопленное произведение вращений Якоби.

Одностороннее вращение в плоскости <math>(j,k)</math> осуществляет процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''':

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

Здесь <math>R(j,k,\theta)</math> - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 244-245)</ref>.

                                                                                
<math>
\begin{matrix}
& & & & j & & k & & & \\
\end{matrix}
</math>

<math>
J_i = 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>a_{jk}</math> не слишком мал"''. Подходящим критерием останова является условие
<math>\begin{align}
|a_{jk}| \ge \varepsilon\sqrt{a_{jj}a_{kk}}
\end{align}</math>.

=== Вычислительное ядро алгоритма ===

Вычислительное ядро составляют множественные вычисления элементов матрицы <math>a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk}</math> и <math>a_{kk} = (G^TG)_{kk}</math>, а также вычисление промежуточной матрицы <math>GR(j, k, \theta)</math> путем умножения исходной матрицы на матрицу вращения Якоби в процессе выполнения процедуры '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''.

После приведения матрицы <math>G^TG</math> к виду, достаточно близкому к диагональному, вычисляются сингулярные числа:

<math>\sigma_i = ||G(:,i)||_2</math>

Матрица левых сингулярных векторов <math>U = [u_1,..., u_n]</math>:

<math>u_i = G(:,i)/\sigma_i</math>

И (если требуется) матрица правых сингулярных векторов <math>V</math> путем накопления произведения вращений Якоби:

<math>
V = J,
J = JR(j, k, \theta)
</math>

=== Макроструктура алгоритма ===
Основную часть метода составляет процедура односторонних вращений матрицы (процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''). Рассмотрим подробнее структуру этой процедуры.

- Вычисление <math>a_{jk}, a_{jj}, a_{kk} (a_{ij} = (G^TG)_{ij})</math>. В связи с тем, что нет необходимости каждый раз вычислять произведение <math>G^TG</math>, для вычисления каждого элемента требуется <math>n</math> операций умножения и <math>(n-1)</math> операций сложения.

- Вычисление <math>\tau, c, s, t</math> не представляют значительной вычислительной сложности.

- Одностороннее вращение матрицы (умножение входной матрицы на матрицу вращения) требует порядка <math>O(n)</math> операций.

В итоге, процедура '''''One-Sided-Jacobi-Rotation''''' требует <math>O(n)</math> операций.

После выполнения алгоритма односторонних вращений необходимо вычислить искомые матрицы <math>S, U, V</math> (матрица сингулярных чисел и матрицы левых и правых векторов).

- Вычисление сингулярных чисел <math>\sigma_i</math> требует <math>n</math> операций умножения.

- Вычисление <math>U</math> требует <math>O(n^2)</math> операций. Аналогично для матрицы <math>V</math>.

=== Схема реализации последовательного алгоритма ===
Метод можно описать следующим образом<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-263)</ref>:

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> - накопленное произведение вращений Якоби.

Метод Якоби для нахождения сингулярных чисел и векторов матрицы <math>G</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

=== Последовательная сложность алгоритма ===
Для вычисления сингулярных чисел и векторов матрицы порядка <math>n</math> в последовательном варианте на каждой итерации алгоритма вызывается процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''', сложность которой равна <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>j</math> — меняется в диапазоне от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения;
* <math>k</math> — меняется в диапазоне от <math>j+1</math> до <math>n</math>, принимая все целочисленные значения.

'''Вторая''' группа вершин соответствует вычислению значений <math>c</math> и <math>s</math> для матрицы вращения

'''Третья''' группа вершин соответствует вычислению произведений матриц <math>GR(k,j, \theta)</math> и <math>JR(k,j, \theta)</math>

<center>
<div class="thumb">
<div class="thumbinner" style="width:{{#expr: 3 * (400 + 35) + 4 * (3 - 1) + 8}}px">
<gallery widths=400px heights=400px>
File:JacobiSVG.png|Рисунок 1. Граф одной итерации алгоритма без отображения входных и выходных данных.
</gallery>
<gallery widths=400px heights=1200px>
File:JacobiSVD2.png|Рисунок 2. Граф одной итерации алгоритма с вычислением матрицы вращения <math>R</math> отображением входных и выходных данных.
</gallery>
</div>
</div>
</center>

=== Ресурс параллелизма алгоритма ===
В данном алгоритме существует возможность параллельного вычисления вращений Якоби<ref>Halil Snopce, Ilir Spahiu Parallel Computation of the SVD, pp. 458-460, http://www.wseas.us/e-library/conferences/2011/Paris/ECC/ECC-73.pdf</ref><ref>V. Strumpen, H. Hoffmann, A. Agarwal A Stream Algorithm for the SVD, pp. 6-7, 23-24, http://publications.csail.mit.edu/lcs/pubs/pdf/MIT-LCS-TM-641.pdf</ref>. Максимальное количество вращений в данном методе равно <math>\frac{n(n-1)}{2}</math>. В параллельной реализации, если <math>n</math> - четное, можно выполнить <math>\frac{n}{2}</math> вращений параллельно, если <math>n</math> - нечетное, то <math>\frac{n-1}{2}</math> вращений. Также следует учесть, что сложность самой процедуры вращения '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''' равна <math>O(n)</math>.

Таким образом, при классификации по высоте ЯПФ метод Якоби относится к алгоритмам со сложностью <math>O(n^2)</math>, а при классификации по ширине ЯПФ — к алгоритмам со сложностью <math>O(n)</math>.

=== Входные и выходные данные алгоритма ===
'''Входные данные''': плотная матрица <math>G</math> (элементы <math>g_{ij}</math>).
Дополнительные ограничения:
* <math>A= G^TG</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

'''Объём входных данных''': <math>n^2</math>

'''Выходные данные''': матрица <math>\Sigma = diag(\sigma_{i})</math>, где <math>\sigma_{i}</math> - сингулярные числа, матрица <math>U</math> левых сингулярных векторов и матрица <math>V</math> правых сингулярных векторов.

'''Объём выходных данных''': <math>2n^2 + n</math> (так как матрица сингулярных чисел <math>\Sigma</math> - диагональная, то достаточно хранить вектор сингулярных чисел (элементы <math>\sigma_i</math>), а также две матрицы <math>U, V</math> левых и правых сингулярных векторов).

=== Свойства алгоритма ===
Метод Якоби является самым медленным из имеющихся алгоритмов вычисления сингулярных чисел и сингулярных векторов матрицы. Тем не менее, для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. Кроме того, в методе Якоби заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, за счет этого метод достигает высокой относительной точности.

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов для метода Якоби является линейным, вычислительная мощность алгоритма также линейна.

Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности (число итераций зависит от входных данных и порогового значения).

== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===

=== Локальность данных и вычислений ===

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

=== Масштабируемость алгоритма и его реализации ===
Данный раздел будет дополнен позже в связи с продлением срока сдачи (до 15 ноября).

=== Динамические характеристики и эффективность реализации алгоритма ===

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

=== Существующие реализации алгоритма ===
Несмотря на то, что метод Якоби является самым медленным из методов вычисления сингулярных чисел и векторов матрицы, он реализован в некоторых пакетах, например, LAPACK, ScaLAPACK, Intel MKL<ref>Intel® Math Kernel Library, https://software.intel.com/ru-ru/node/521153#CAF8CA08-1C69-48F8-B0C3-B2B27AFA92FA</ref>.

== Литература ==

Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов

2016-10-27T14:18:19Z

Светлана: /* Информационный граф */

{{algorithm
| name = Метод Якоби вычисления сингулярных чисел и векторов
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math>O(n^2)</math>
| pf_width = <math>O(n)</math>
| input_data = <math>n^2</math>
| output_data = <math>2n^2 + n</math>
}}
Авторы статьи: [[Участник:Почернина Елена|Почернина Елена]] (группа 601), [[Участница:Светлана|Костюкова Светлана]] (группа 601)

== Свойства и структура алгоритма ==

=== Общее описание алгоритма ===

'''Метод Якоби''' позволяет найти сингулярное разложение плотной матрицы <math>G</math> размера <math>n</math>×<math>n</math>. Такое разложение можно записать в виде <math>G = U\Sigma V^*</math>, где <math>U, V</math> - унитарные (''в вещественном случае - ортогональные'') матрицы размера <math>n</math>×<math>n</math>, <math>\Sigma</math> - диагональная матрица размера <math>n</math>×<math>n</math> с вещественными положительными числами на главной диагонали. Диагональные элементы матрицы <math>\Sigma</math> называются сингулярными числами матрицы <math>G</math>, а столбцы матриц <math>U,V</math> - левыми и правыми сингулярными векторами.

''Неявный''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref> метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице <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>. Поэтому алгоритм называется ''методом односторонних вращений''.

Метод Якоби для вычисления сингулярных чисел и векторов является самым медленным из имеющихся алгоритмов поиска сингулярных чисел и векторов, но тем не менее, интерес к нему сохраняется. Для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. В частности, метод Якоби вычисляет сингулярные числа матрицы <math>G</math> с высокой точностью, если G может быть представлена в виде <math>G = DX</math>, где <math>D</math> – диагональная матрица, а <math>X</math> – хорошо обусловлена. В этом случае заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, в то время как другие алгоритмы включают в себя приведение матрицы к двухдиагональной форме, из-за чего и теряют все верные разряды во всех сингулярных числах, кроме старшего.

=== Математическое описание алгоритма ===
Исходные данные: матрица <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> - матрица правых сингулярных векторов.

Основной алгоритм заключается в следующем<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref>:

Для всех элементов матрицы <math>A^{(i)} = (G^TG)^{(i)}</math>, находящихся вне главной диагонали <math>(a_{12}, a_{13},\dots, a_{1n}, a_{23}, a_{24},\dots, a_{2n},\dots, a_{n-1n})</math> проверяем условие <math>|A_{jk}|<\varepsilon</math>:
Если <math>|A_{jk}^{(i)}| \le \varepsilon</math>, то переходим к следующему элементу,
Если <math>|A_{jk}^{(i)}| > \varepsilon</math>, то нужно выполнить вращение в плоскости <math>(j,k)</math>, после чего перейти к следующей итерации алгоритма.
Если на очередной итерации не было совершено ни одного вращения, значит матрица <math>A</math> стала достаточна близка к диагональной, и можно перейти к вычислению искомых значений:

<math>\sigma_i = ||G(:,i)||_2</math>;

<math>U = [u_1,\dots,u_n]</math>, где <math>u_i = G(:,i)/\sigma_i</math>;

<math>V = J</math>, где <math>J</math> - накопленное произведение вращений Якоби.

Одностороннее вращение в плоскости <math>(j,k)</math> осуществляет процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''':

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

Здесь <math>R(j,k,\theta)</math> - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 244-245)</ref>.

                                                                                
<math>
\begin{matrix}
& & & & j & & k & & & \\
\end{matrix}
</math>

<math>
J_i = 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>a_{jk}</math> не слишком мал"''. Подходящим критерием останова является условие
<math>\begin{align}
|a_{jk}| \ge \varepsilon\sqrt{a_{jj}a_{kk}}
\end{align}</math>.

=== Вычислительное ядро алгоритма ===

Вычислительное ядро составляют множественные вычисления элементов матрицы <math>a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk}</math> и <math>a_{kk} = (G^TG)_{kk}</math>, а также вычисление промежуточной матрицы <math>GR(j, k, \theta)</math> путем умножения исходной матрицы на матрицу вращения Якоби в процессе выполнения процедуры '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''.

После приведения матрицы <math>G^TG</math> к виду, достаточно близкому к диагональному, вычисляются сингулярные числа:

<math>\sigma_i = ||G(:,i)||_2</math>

Матрица левых сингулярных векторов <math>U = [u_1,..., u_n]</math>:

<math>u_i = G(:,i)/\sigma_i</math>

И (если требуется) матрица правых сингулярных векторов <math>V</math> путем накопления произведения вращений Якоби:

<math>
V = J,
J = JR(j, k, \theta)
</math>

=== Макроструктура алгоритма ===
Основную часть метода составляет процедура односторонних вращений матрицы (процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''). Рассмотрим подробнее структуру этой процедуры.

- Вычисление <math>a_{jk}, a_{jj}, a_{kk} (a_{ij} = (G^TG)_{ij})</math>. В связи с тем, что нет необходимости каждый раз вычислять произведение <math>G^TG</math>, для вычисления каждого элемента требуется <math>n</math> операций умножения и <math>(n-1)</math> операций сложения.

- Вычисление <math>\tau, c, s, t</math> не представляют значительной вычислительной сложности.

- Одностороннее вращение матрицы (умножение входной матрицы на матрицу вращения) требует порядка <math>O(n)</math> операций.

В итоге, процедура '''''One-Sided-Jacobi-Rotation''''' требует <math>O(n)</math> операций.

После выполнения алгоритма односторонних вращений необходимо вычислить искомые матрицы <math>S, U, V</math> (матрица сингулярных чисел и матрицы левых и правых векторов).

- Вычисление сингулярных чисел <math>\sigma_i</math> требует <math>n</math> операций умножения.

- Вычисление <math>U</math> требует <math>O(n^2)</math> операций. Аналогично для матрицы <math>V</math>.

=== Схема реализации последовательного алгоритма ===
Метод можно описать следующим образом<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-263)</ref>:

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> - накопленное произведение вращений Якоби.

Метод Якоби для нахождения сингулярных чисел и векторов матрицы <math>G</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

=== Последовательная сложность алгоритма ===
Для вычисления сингулярных чисел и векторов матрицы порядка <math>n</math> в последовательном варианте на каждой итерации алгоритма вызывается процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''', сложность которой равна <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>j</math> — меняется в диапазоне от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения;
* <math>k</math> — меняется в диапазоне от <math>j+1</math> до <math>n</math>, принимая все целочисленные значения.

'''Вторая''' группа вершин соответствует вычислению значений <math>c</math> и <math>s</math> для матрицы вращения

'''Третья''' группа вершин соответствует вычислению произведений матриц <math>GR(k,j, \theta)</math> и <math>JR(k,j, \theta)</math>

<center>
<div class="thumb">
<div class="thumbinner" style="width:{{#expr: 3 * (400 + 35) + 4 * (3 - 1) + 8}}px">
<gallery widths=400px heights=400px>
File:JacobiSVG.png|Рисунок 1. Граф одной итерации алгоритма без отображения входных и выходных данных.
</gallery>
<gallery widths=400px heights=1200px>
File:JacobiSVD2.png|Рисунок 2. Граф одной итерации алгоритма с отображением входных и выходных данных.
</gallery>
</div>
</div>
</center>

=== Ресурс параллелизма алгоритма ===
В данном алгоритме существует возможность параллельного вычисления вращений Якоби<ref>Halil Snopce, Ilir Spahiu Parallel Computation of the SVD, pp. 458-460, http://www.wseas.us/e-library/conferences/2011/Paris/ECC/ECC-73.pdf</ref><ref>V. Strumpen, H. Hoffmann, A. Agarwal A Stream Algorithm for the SVD, pp. 6-7, 23-24, http://publications.csail.mit.edu/lcs/pubs/pdf/MIT-LCS-TM-641.pdf</ref>. Максимальное количество вращений в данном методе равно <math>\frac{n(n-1)}{2}</math>. В параллельной реализации, если <math>n</math> - четное, можно выполнить <math>\frac{n}{2}</math> вращений параллельно, если <math>n</math> - нечетное, то <math>\frac{n-1}{2}</math> вращений. Также следует учесть, что сложность самой процедуры вращения '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''' равна <math>O(n)</math>.

Таким образом, при классификации по высоте ЯПФ метод Якоби относится к алгоритмам со сложностью <math>O(n^2)</math>, а при классификации по ширине ЯПФ — к алгоритмам со сложностью <math>O(n)</math>.

=== Входные и выходные данные алгоритма ===
'''Входные данные''': плотная матрица <math>G</math> (элементы <math>g_{ij}</math>).
Дополнительные ограничения:
* <math>A= G^TG</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

'''Объём входных данных''': <math>n^2</math>

'''Выходные данные''': матрица <math>\Sigma = diag(\sigma_{i})</math>, где <math>\sigma_{i}</math> - сингулярные числа, матрица <math>U</math> левых сингулярных векторов и матрица <math>V</math> правых сингулярных векторов.

'''Объём выходных данных''': <math>2n^2 + n</math> (так как матрица сингулярных чисел <math>\Sigma</math> - диагональная, то достаточно хранить вектор сингулярных чисел (элементы <math>\sigma_i</math>), а также две матрицы <math>U, V</math> левых и правых сингулярных векторов).

=== Свойства алгоритма ===
Метод Якоби является самым медленным из имеющихся алгоритмов вычисления сингулярных чисел и сингулярных векторов матрицы. Тем не менее, для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. Кроме того, в методе Якоби заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, за счет этого метод достигает высокой относительной точности.

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов для метода Якоби является линейным, вычислительная мощность алгоритма также линейна.

Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности (число итераций зависит от входных данных и порогового значения).

== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===

=== Локальность данных и вычислений ===

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

=== Масштабируемость алгоритма и его реализации ===
Данный раздел будет дополнен позже в связи с продлением срока сдачи (до 15 ноября).

=== Динамические характеристики и эффективность реализации алгоритма ===

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

=== Существующие реализации алгоритма ===
Несмотря на то, что метод Якоби является самым медленным из методов вычисления сингулярных чисел и векторов матрицы, он реализован в некоторых пакетах, например, LAPACK, ScaLAPACK, Intel MKL<ref>Intel® Math Kernel Library, https://software.intel.com/ru-ru/node/521153#CAF8CA08-1C69-48F8-B0C3-B2B27AFA92FA</ref>.

== Литература ==

Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов

2016-10-27T14:16:44Z

Светлана: /* Информационный граф */

{{algorithm
| name = Метод Якоби вычисления сингулярных чисел и векторов
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math>O(n^2)</math>
| pf_width = <math>O(n)</math>
| input_data = <math>n^2</math>
| output_data = <math>2n^2 + n</math>
}}
Авторы статьи: [[Участник:Почернина Елена|Почернина Елена]] (группа 601), [[Участница:Светлана|Костюкова Светлана]] (группа 601)

== Свойства и структура алгоритма ==

=== Общее описание алгоритма ===

'''Метод Якоби''' позволяет найти сингулярное разложение плотной матрицы <math>G</math> размера <math>n</math>×<math>n</math>. Такое разложение можно записать в виде <math>G = U\Sigma V^*</math>, где <math>U, V</math> - унитарные (''в вещественном случае - ортогональные'') матрицы размера <math>n</math>×<math>n</math>, <math>\Sigma</math> - диагональная матрица размера <math>n</math>×<math>n</math> с вещественными положительными числами на главной диагонали. Диагональные элементы матрицы <math>\Sigma</math> называются сингулярными числами матрицы <math>G</math>, а столбцы матриц <math>U,V</math> - левыми и правыми сингулярными векторами.

''Неявный''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref> метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице <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>. Поэтому алгоритм называется ''методом односторонних вращений''.

Метод Якоби для вычисления сингулярных чисел и векторов является самым медленным из имеющихся алгоритмов поиска сингулярных чисел и векторов, но тем не менее, интерес к нему сохраняется. Для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. В частности, метод Якоби вычисляет сингулярные числа матрицы <math>G</math> с высокой точностью, если G может быть представлена в виде <math>G = DX</math>, где <math>D</math> – диагональная матрица, а <math>X</math> – хорошо обусловлена. В этом случае заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, в то время как другие алгоритмы включают в себя приведение матрицы к двухдиагональной форме, из-за чего и теряют все верные разряды во всех сингулярных числах, кроме старшего.

=== Математическое описание алгоритма ===
Исходные данные: матрица <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> - матрица правых сингулярных векторов.

Основной алгоритм заключается в следующем<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref>:

Для всех элементов матрицы <math>A^{(i)} = (G^TG)^{(i)}</math>, находящихся вне главной диагонали <math>(a_{12}, a_{13},\dots, a_{1n}, a_{23}, a_{24},\dots, a_{2n},\dots, a_{n-1n})</math> проверяем условие <math>|A_{jk}|<\varepsilon</math>:
Если <math>|A_{jk}^{(i)}| \le \varepsilon</math>, то переходим к следующему элементу,
Если <math>|A_{jk}^{(i)}| > \varepsilon</math>, то нужно выполнить вращение в плоскости <math>(j,k)</math>, после чего перейти к следующей итерации алгоритма.
Если на очередной итерации не было совершено ни одного вращения, значит матрица <math>A</math> стала достаточна близка к диагональной, и можно перейти к вычислению искомых значений:

<math>\sigma_i = ||G(:,i)||_2</math>;

<math>U = [u_1,\dots,u_n]</math>, где <math>u_i = G(:,i)/\sigma_i</math>;

<math>V = J</math>, где <math>J</math> - накопленное произведение вращений Якоби.

Одностороннее вращение в плоскости <math>(j,k)</math> осуществляет процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''':

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

Здесь <math>R(j,k,\theta)</math> - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 244-245)</ref>.

                                                                                
<math>
\begin{matrix}
& & & & j & & k & & & \\
\end{matrix}
</math>

<math>
J_i = 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>a_{jk}</math> не слишком мал"''. Подходящим критерием останова является условие
<math>\begin{align}
|a_{jk}| \ge \varepsilon\sqrt{a_{jj}a_{kk}}
\end{align}</math>.

=== Вычислительное ядро алгоритма ===

Вычислительное ядро составляют множественные вычисления элементов матрицы <math>a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk}</math> и <math>a_{kk} = (G^TG)_{kk}</math>, а также вычисление промежуточной матрицы <math>GR(j, k, \theta)</math> путем умножения исходной матрицы на матрицу вращения Якоби в процессе выполнения процедуры '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''.

После приведения матрицы <math>G^TG</math> к виду, достаточно близкому к диагональному, вычисляются сингулярные числа:

<math>\sigma_i = ||G(:,i)||_2</math>

Матрица левых сингулярных векторов <math>U = [u_1,..., u_n]</math>:

<math>u_i = G(:,i)/\sigma_i</math>

И (если требуется) матрица правых сингулярных векторов <math>V</math> путем накопления произведения вращений Якоби:

<math>
V = J,
J = JR(j, k, \theta)
</math>

=== Макроструктура алгоритма ===
Основную часть метода составляет процедура односторонних вращений матрицы (процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''). Рассмотрим подробнее структуру этой процедуры.

- Вычисление <math>a_{jk}, a_{jj}, a_{kk} (a_{ij} = (G^TG)_{ij})</math>. В связи с тем, что нет необходимости каждый раз вычислять произведение <math>G^TG</math>, для вычисления каждого элемента требуется <math>n</math> операций умножения и <math>(n-1)</math> операций сложения.

- Вычисление <math>\tau, c, s, t</math> не представляют значительной вычислительной сложности.

- Одностороннее вращение матрицы (умножение входной матрицы на матрицу вращения) требует порядка <math>O(n)</math> операций.

В итоге, процедура '''''One-Sided-Jacobi-Rotation''''' требует <math>O(n)</math> операций.

После выполнения алгоритма односторонних вращений необходимо вычислить искомые матрицы <math>S, U, V</math> (матрица сингулярных чисел и матрицы левых и правых векторов).

- Вычисление сингулярных чисел <math>\sigma_i</math> требует <math>n</math> операций умножения.

- Вычисление <math>U</math> требует <math>O(n^2)</math> операций. Аналогично для матрицы <math>V</math>.

=== Схема реализации последовательного алгоритма ===
Метод можно описать следующим образом<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-263)</ref>:

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> - накопленное произведение вращений Якоби.

Метод Якоби для нахождения сингулярных чисел и векторов матрицы <math>G</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

=== Последовательная сложность алгоритма ===
Для вычисления сингулярных чисел и векторов матрицы порядка <math>n</math> в последовательном варианте на каждой итерации алгоритма вызывается процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''', сложность которой равна <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>j</math> — меняется в диапазоне от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения;
* <math>k</math> — меняется в диапазоне от <math>j+1</math> до <math>n</math>, принимая все целочисленные значения.

'''Вторая''' группа вершин соответствует вычислению значений <math>c</math> и <math>s</math> для матрицы вращения

'''Третья''' группа вершин соответствует вычислению произведений матриц <math>GR(k,j, \theta)</math> и <math>JR(k,j, \theta)</math>

<center>
<div class="thumb">
<div class="thumbinner" style="width:{{#expr: 3 * (400 + 35) + 4 * (3 - 1) + 8}}px">
<gallery widths=400px heights=1200px>
File:JacobiSVG.png|Рисунок 1. Граф одной итерации алгоритма без отображения входных и выходных данных.
File:JacobiSVD2.png|Рисунок 2. Граф одной итерации алгоритма с отображением входных и выходных данных.
</gallery>
</div>
</div>
</center>

=== Ресурс параллелизма алгоритма ===
В данном алгоритме существует возможность параллельного вычисления вращений Якоби<ref>Halil Snopce, Ilir Spahiu Parallel Computation of the SVD, pp. 458-460, http://www.wseas.us/e-library/conferences/2011/Paris/ECC/ECC-73.pdf</ref><ref>V. Strumpen, H. Hoffmann, A. Agarwal A Stream Algorithm for the SVD, pp. 6-7, 23-24, http://publications.csail.mit.edu/lcs/pubs/pdf/MIT-LCS-TM-641.pdf</ref>. Максимальное количество вращений в данном методе равно <math>\frac{n(n-1)}{2}</math>. В параллельной реализации, если <math>n</math> - четное, можно выполнить <math>\frac{n}{2}</math> вращений параллельно, если <math>n</math> - нечетное, то <math>\frac{n-1}{2}</math> вращений. Также следует учесть, что сложность самой процедуры вращения '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''' равна <math>O(n)</math>.

Таким образом, при классификации по высоте ЯПФ метод Якоби относится к алгоритмам со сложностью <math>O(n^2)</math>, а при классификации по ширине ЯПФ — к алгоритмам со сложностью <math>O(n)</math>.

=== Входные и выходные данные алгоритма ===
'''Входные данные''': плотная матрица <math>G</math> (элементы <math>g_{ij}</math>).
Дополнительные ограничения:
* <math>A= G^TG</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

'''Объём входных данных''': <math>n^2</math>

'''Выходные данные''': матрица <math>\Sigma = diag(\sigma_{i})</math>, где <math>\sigma_{i}</math> - сингулярные числа, матрица <math>U</math> левых сингулярных векторов и матрица <math>V</math> правых сингулярных векторов.

'''Объём выходных данных''': <math>2n^2 + n</math> (так как матрица сингулярных чисел <math>\Sigma</math> - диагональная, то достаточно хранить вектор сингулярных чисел (элементы <math>\sigma_i</math>), а также две матрицы <math>U, V</math> левых и правых сингулярных векторов).

=== Свойства алгоритма ===
Метод Якоби является самым медленным из имеющихся алгоритмов вычисления сингулярных чисел и сингулярных векторов матрицы. Тем не менее, для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. Кроме того, в методе Якоби заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, за счет этого метод достигает высокой относительной точности.

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов для метода Якоби является линейным, вычислительная мощность алгоритма также линейна.

Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности (число итераций зависит от входных данных и порогового значения).

== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===

=== Локальность данных и вычислений ===

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

=== Масштабируемость алгоритма и его реализации ===
Данный раздел будет дополнен позже в связи с продлением срока сдачи (до 15 ноября).

=== Динамические характеристики и эффективность реализации алгоритма ===

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

=== Существующие реализации алгоритма ===
Несмотря на то, что метод Якоби является самым медленным из методов вычисления сингулярных чисел и векторов матрицы, он реализован в некоторых пакетах, например, LAPACK, ScaLAPACK, Intel MKL<ref>Intel® Math Kernel Library, https://software.intel.com/ru-ru/node/521153#CAF8CA08-1C69-48F8-B0C3-B2B27AFA92FA</ref>.

== Литература ==

Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов

2016-10-15T11:07:06Z

Светлана: /* Информационный граф */

{{algorithm
| name = Метод Якоби вычисления сингулярных чисел и векторов
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math>O(n^2)</math>
| pf_width = <math>O(n)</math>
| input_data = <math>n^2</math>
| output_data = <math>2n^2 + n</math>
}}
Авторы статьи: [[Участник:Почернина Елена|Почернина Елена]] (группа 601), [[Участница:Светлана|Костюкова Светлана]] (группа 601)

== Свойства и структура алгоритма ==

=== Общее описание алгоритма ===

'''Метод Якоби''' позволяет найти сингулярное разложение плотной матрицы <math>G</math> размера <math>n</math>×<math>n</math>. Такое разложение можно записать в виде <math>G = U\Sigma V^*</math>, где <math>U, V</math> - унитарные (''в вещественном случае - ортогональные'') матрицы размера <math>n</math>×<math>n</math>, <math>\Sigma</math> - диагональная матрица размера <math>n</math>×<math>n</math> с вещественными положительными числами на главной диагонали. Диагональные элементы матрицы <math>\Sigma</math> называются сингулярными числами матрицы <math>G</math>, а столбцы матриц <math>U,V</math> - левыми и правыми сингулярными векторами.

''Неявный''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref> метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице <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>. Поэтому алгоритм называется ''методом односторонних вращений''.

Метод Якоби для вычисления сингулярных чисел и векторов является самым медленным из имеющихся алгоритмов поиска сингулярных чисел и векторов, но тем не менее, интерес к нему сохраняется. Для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. В частности, метод Якоби вычисляет сингулярные числа матрицы <math>G</math> с высокой точностью, если G может быть представлена в виде <math>G = DX</math>, где <math>D</math> – диагональная матрица, а <math>X</math> – хорошо обусловлена. В этом случае заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, в то время как другие алгоритмы включают в себя приведение матрицы к двухдиагональной форме, из-за чего и теряют все верные разряды во всех сингулярных числах, кроме старшего.

=== Математическое описание алгоритма ===
Исходные данные: матрица <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> - матрица правых сингулярных векторов.

Основной алгоритм заключается в следующем<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref>:

Для всех элементов матрицы <math>A^{(i)} = (G^TG)^{(i)}</math>, находящихся вне главной диагонали <math>(a_{12}, a_{13},\dots, a_{1n}, a_{23}, a_{24},\dots, a_{2n},\dots, a_{n-1n})</math> проверяем условие <math>|A_{jk}|<\varepsilon</math>:
Если <math>|A_{jk}^{(i)}| \le \varepsilon</math>, то переходим к следующему элементу,
Если <math>|A_{jk}^{(i)}| > \varepsilon</math>, то нужно выполнить вращение в плоскости <math>(j,k)</math>, после чего перейти к следующей итерации алгоритма.
Если на очередной итерации не было совершено ни одного вращения, значит матрица <math>A</math> стала достаточна близка к диагональной, и можно перейти к вычислению искомых значений:

<math>\sigma_i = ||G(:,i)||_2</math>;

<math>U = [u_1,\dots,u_n]</math>, где <math>u_i = G(:,i)/\sigma_i</math>;

<math>V = J</math>, где <math>J</math> - накопленное произведение вращений Якоби.

Одностороннее вращение в плоскости <math>(j,k)</math> осуществляет процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''':

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

Здесь <math>R(j,k,\theta)</math> - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 244-245)</ref>.

                                                                                
<math>
\begin{matrix}
& & & & j & & k & & & \\
\end{matrix}
</math>

<math>
J_i = 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>a_{jk}</math> не слишком мал"''. Подходящим критерием останова является условие
<math>\begin{align}
|a_{jk}| \ge \varepsilon\sqrt{a_{jj}a_{kk}}
\end{align}</math>.

=== Вычислительное ядро алгоритма ===

Вычислительное ядро составляют множественные вычисления элементов матрицы <math>a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk}</math> и <math>a_{kk} = (G^TG)_{kk}</math>, а также вычисление промежуточной матрицы <math>GR(j, k, \theta)</math> путем умножения исходной матрицы на матрицу вращения Якоби в процессе выполнения процедуры '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''.

После приведения матрицы <math>G^TG</math> к виду, достаточно близкому к диагональному, вычисляются сингулярные числа:

<math>\sigma_i = ||G(:,i)||_2</math>

Матрица левых сингулярных векторов <math>U = [u_1,..., u_n]</math>:

<math>u_i = G(:,i)/\sigma_i</math>

И (если требуется) матрица правых сингулярных векторов <math>V</math> путем накопления произведения вращений Якоби:

<math>
V = J,
J = JR(j, k, \theta)
</math>

=== Макроструктура алгоритма ===
Основную часть метода составляет процедура односторонних вращений матрицы (процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''). Рассмотрим подробнее структуру этой процедуры.

- Вычисление <math>a_{jk}, a_{jj}, a_{kk} (a_{ij} = (G^TG)_{ij})</math>. В связи с тем, что нет необходимости каждый раз вычислять произведение <math>G^TG</math>, для вычисления каждого элемента требуется <math>n</math> операций умножения и <math>(n-1)</math> операций сложения.

- Вычисление <math>\tau, c, s, t</math> не представляют значительной вычислительной сложности.

- Одностороннее вращение матрицы (умножение входной матрицы на матрицу вращения) требует порядка <math>O(n)</math> операций.

В итоге, процедура '''''One-Sided-Jacobi-Rotation''''' требует <math>O(n)</math> операций.

После выполнения алгоритма односторонних вращений необходимо вычислить искомые матрицы <math>S, U, V</math> (матрица сингулярных чисел и матрицы левых и правых векторов).

- Вычисление сингулярных чисел <math>\sigma_i</math> требует <math>n</math> операций умножения.

- Вычисление <math>U</math> требует <math>O(n^2)</math> операций. Аналогично для матрицы <math>V</math>.

=== Схема реализации последовательного алгоритма ===
Метод можно описать следующим образом<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-263)</ref>:

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> - накопленное произведение вращений Якоби.

Метод Якоби для нахождения сингулярных чисел и векторов матрицы <math>G</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

=== Последовательная сложность алгоритма ===
Для вычисления сингулярных чисел и векторов матрицы порядка <math>n</math> в последовательном варианте на каждой итерации алгоритма вызывается процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''', сложность которой равна <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>j</math> — меняется в диапазоне от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения;
* <math>k</math> — меняется в диапазоне от <math>j+1</math> до <math>n</math>, принимая все целочисленные значения.

'''Вторая''' группа вершин соответствует вычислению значений <math>c</math> и <math>s</math> для матрицы вращения

'''Третья''' группа вершин соответствует вычислению произведений матриц <math>GR(k,j, \theta)</math> и <math>JR(k,j, \theta)</math>

<center>
<div class="thumb">
<div class="thumbinner" style="width:{{#expr: 3 * (400 + 35) + 4 * (3 - 1) + 8}}px">
<gallery widths=400px heights=400px>
File:JacobiSVG.png|Рисунок 1. Граф одной итерации алгоритма без отображения входных и выходных данных.
File:JacobiSVD2.png|Рисунок 2. Граф одной итерации алгоритма с отображением входных и выходных данных.
</gallery>
</div>
</div>
</center>

=== Ресурс параллелизма алгоритма ===
В данном алгоритме существует возможность параллельного вычисления вращений Якоби<ref>Halil Snopce, Ilir Spahiu Parallel Computation of the SVD, pp. 458-460, http://www.wseas.us/e-library/conferences/2011/Paris/ECC/ECC-73.pdf</ref>. Максимальное количество вращений в данном методе равно <math>\frac{n(n-1)}{2}</math>. В параллельной реализации, если <math>n</math> - четное, можно выполнить <math>\frac{n}{2}</math> вращений параллельно, если <math>n</math> - нечетное, то <math>\frac{n-1}{2}</math> вращений. Также следует учесть, что сложность самой процедуры вращения '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''' равна <math>O(n)</math>.

Таким образом, при классификации по высоте ЯПФ метод Якоби относится к алгоритмам со сложностью <math>O(n^2)</math>, а при классификации по ширине ЯПФ — к алгоритмам со сложностью <math>O(n)</math>.

=== Входные и выходные данные алгоритма ===
'''Входные данные''': плотная матрица <math>G</math> (элементы <math>g_{ij}</math>).
Дополнительные ограничения:
* <math>A= G^TG</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

'''Объём входных данных''': <math>n^2</math>

'''Выходные данные''': матрица <math>\Sigma = diag(\sigma_{i})</math>, где <math>\sigma_{i}</math> - сингулярные числа, матрица <math>U</math> левых сингулярных векторов и матрица <math>V</math> правых сингулярных векторов.

'''Объём выходных данных''': <math>2n^2 + n</math> (так как матрица сингулярных чисел <math>\Sigma</math> - диагональная, то достаточно хранить вектор сингулярных чисел (элементы <math>\sigma_i</math>), а также две матрицы <math>U, V</math> левых и правых сингулярных векторов).

=== Свойства алгоритма ===
Метод Якоби является самым медленным из имеющихся алгоритмов вычисления сингулярных чисел и сингулярных векторов матрицы. Тем не менее, для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. Кроме того, в методе Якоби заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, за счет этого метод достигает высокой относительной точности.

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов для метода Якоби является линейным, вычислительная мощность алгоритма также линейна.

Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности (число итераций зависит от входных данных и порогового значения).

== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===

=== Локальность данных и вычислений ===

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

=== Масштабируемость алгоритма и его реализации ===
Данный раздел будет дополнен позже в связи с продлением срока сдачи (до 15 ноября).

=== Динамические характеристики и эффективность реализации алгоритма ===

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

=== Существующие реализации алгоритма ===
Несмотря на то, что метод Якоби является самым медленным из методов вычисления сингулярных чисел и векторов матрицы, он реализован в некоторых пакетах, например, LAPACK,scaLAPACK, Intel MKL<ref>Intel® Math Kernel Library, https://software.intel.com/ru-ru/node/521153#CAF8CA08-1C69-48F8-B0C3-B2B27AFA92FA</ref>.

== Литература ==

Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов

2016-10-15T10:12:05Z

Светлана: /* Информационный граф */

{{algorithm
| name = Метод Якоби вычисления сингулярных чисел и векторов
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math>O(n^2)</math>
| pf_width = <math>O(n)</math>
| input_data = <math>n^2</math>
| output_data = <math>2n^2 + n</math>
}}
Авторы статьи: [[Участник:Почернина Елена|Почернина Елена]] (группа 601), [[Участница:Светлана|Костюкова Светлана]] (группа 601)

== Свойства и структура алгоритма ==

=== Общее описание алгоритма ===

'''Метод Якоби''' позволяет найти сингулярное разложение плотной матрицы <math>G</math> размера <math>n</math>×<math>n</math>. Такое разложение можно записать в виде <math>G = U\Sigma V^*</math>, где <math>U, V</math> - унитарные (''в вещественном случае - ортогональные'') матрицы размера <math>n</math>×<math>n</math>, <math>\Sigma</math> - диагональная матрица размера <math>n</math>×<math>n</math> с вещественными положительными числами на главной диагонали. Диагональные элементы матрицы <math>\Sigma</math> называются сингулярными числами матрицы <math>G</math>, а столбцы матриц <math>U,V</math> - левыми и правыми сингулярными векторами.

''Неявный''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref> метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице <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>. Поэтому алгоритм называется ''методом односторонних вращений''.

Метод Якоби для вычисления сингулярных чисел и векторов является самым медленным из имеющихся алгоритмов поиска сингулярных чисел и векторов, но тем не менее, интерес к нему сохраняется. Для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. В частности, метод Якоби вычисляет сингулярные числа матрицы <math>G</math> с высокой точностью, если G может быть представлена в виде <math>G = DX</math>, где <math>D</math> – диагональная матрица, а <math>X</math> – хорошо обусловлена. В этом случае заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, в то время как другие алгоритмы включают в себя приведение матрицы к двухдиагональной форме, из-за чего и теряют все верные разряды во всех сингулярных числах, кроме старшего.

=== Математическое описание алгоритма ===
Исходные данные: матрица <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> - матрица правых сингулярных векторов.

Основной алгоритм заключается в следующем<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref>:

Для всех элементов матрицы <math>A^{(i)} = (G^TG)^{(i)}</math>, находящихся вне главной диагонали <math>(a_{12}, a_{13},\dots, a_{1n}, a_{23}, a_{24},\dots, a_{2n},\dots, a_{n-1n})</math> проверяем условие <math>|A_{jk}|<\varepsilon</math>:
Если <math>|A_{jk}^{(i)}| \le \varepsilon</math>, то переходим к следующему элементу,
Если <math>|A_{jk}^{(i)}| > \varepsilon</math>, то нужно выполнить вращение в плоскости <math>(j,k)</math>, после чего перейти к следующей итерации алгоритма.
Если на очередной итерации не было совершено ни одного вращения, значит матрица <math>A</math> стала достаточна близка к диагональной, и можно перейти к вычислению искомых значений:

<math>\sigma_i = ||G(:,i)||_2</math>;

<math>U = [u_1,\dots,u_n]</math>, где <math>u_i = G(:,i)/\sigma_i</math>;

<math>V = J</math>, где <math>J</math> - накопленное произведение вращений Якоби.

Одностороннее вращение в плоскости <math>(j,k)</math> осуществляет процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''':

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

Здесь <math>R(j,k,\theta)</math> - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 244-245)</ref>.

                                                                                
<math>
\begin{matrix}
& & & & j & & k & & & \\
\end{matrix}
</math>

<math>
J_i = 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>a_{jk}</math> не слишком мал"''. Подходящим критерием останова является условие
<math>\begin{align}
|a_{jk}| \ge \varepsilon\sqrt{a_{jj}a_{kk}}
\end{align}</math>.

=== Вычислительное ядро алгоритма ===

Вычислительное ядро составляют множественные вычисления элементов матрицы <math>a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk}</math> и <math>a_{kk} = (G^TG)_{kk}</math>, а также вычисление промежуточной матрицы <math>GR(j, k, \theta)</math> путем умножения исходной матрицы на матрицу вращения Якоби в процессе выполнения процедуры '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''.

После приведения матрицы <math>G^TG</math> к виду, достаточно близкому к диагональному, вычисляются сингулярные числа:

<math>\sigma_i = ||G(:,i)||_2</math>

Матрица левых сингулярных векторов <math>U = [u_1,..., u_n]</math>:

<math>u_i = G(:,i)/\sigma_i</math>

И (если требуется) матрица правых сингулярных векторов <math>V</math> путем накопления произведения вращений Якоби:

<math>
V = J,
J = JR(j, k, \theta)
</math>

=== Макроструктура алгоритма ===
Основную часть метода составляет процедура односторонних вращений матрицы (процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''). Рассмотрим подробнее структуру этой процедуры.

- Вычисление <math>a_{jk}, a_{jj}, a_{kk} (a_{ij} = (G^TG)_{ij})</math>. В связи с тем, что нет необходимости каждый раз вычислять произведение <math>G^TG</math>, для вычисления каждого элемента требуется <math>n</math> операций умножения и <math>(n-1)</math> операций сложения.

- Вычисление <math>\tau, c, s, t</math> не представляют значительной вычислительной сложности.

- Одностороннее вращение матрицы (умножение входной матрицы на матрицу вращения) требует порядка <math>O(n)</math> операций.

В итоге, процедура '''''One-Sided-Jacobi-Rotation''''' требует <math>O(n)</math> операций.

После выполнения алгоритма односторонних вращений необходимо вычислить искомые матрицы <math>S, U, V</math> (матрица сингулярных чисел и матрицы левых и правых векторов).

- Вычисление сингулярных чисел <math>\sigma_i</math> требует <math>n</math> операций умножения.

- Вычисление <math>U</math> требует <math>O(n^2)</math> операций. Аналогично для матрицы <math>V</math>.

=== Схема реализации последовательного алгоритма ===
Метод можно описать следующим образом<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-263)</ref>:

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> - накопленное произведение вращений Якоби.

Метод Якоби для нахождения сингулярных чисел и векторов матрицы <math>G</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

=== Последовательная сложность алгоритма ===
Для вычисления сингулярных чисел и векторов матрицы порядка <math>n</math> в последовательном варианте на каждой итерации алгоритма вызывается процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''', сложность которой равна <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>j</math> — меняется в диапазоне от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения;
* <math>k</math> — меняется в диапазоне от <math>j+1</math> до <math>n</math>, принимая все целочисленные значения.

'''Вторая''' группа вершин соответствует вычислению значений <math>c</math> и <math>s</math> для матрицы вращения

'''Третья''' группа вершин соответствует вычислению произведений матриц <math>GR(k,j, \theta)</math> и <math>JR(k,j, \theta)</math>

<center>
<div class="thumb">
<div class="thumbinner" style="width:{{#expr: 3 * (400 + 35) + 4 * (3 - 1) + 8}}px">
<gallery widths=400px heights=400px>
File:JacobiSVG.png|Рисунок 1. Граф одной итерации алгоритма без отображения входных и выходных данных.
File:JacobiSVD2.png|Рисунок 2. Граф одной итерации алгоритма с отображением входных и выходных данных.
</gallery>
</div>
</div>
</center>

=== Ресурс параллелизма алгоритма ===

Таким образом, при классификации по высоте ЯПФ метод Якоби относится к алгоритмам со сложностью <math>O(n^2)</math>, а при классификации по ширине ЯПФ — к алгоритмам со сложностью <math>O(n)</math>.

=== Входные и выходные данные алгоритма ===
'''Входные данные''': плотная матрица <math>G</math> (элементы <math>g_{ij}</math>).
Дополнительные ограничения:
* <math>A= G^TG</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

'''Объём входных данных''': <math>n^2</math>

'''Выходные данные''': матрица <math>\Sigma = diag(\sigma_{i})</math>, где <math>\sigma_{i}</math> - сингулярные числа, матрица <math>U</math> левых сингулярных векторов и матрица <math>V</math> правых сингулярных векторов.

'''Объём выходных данных''': <math>2n^2 + n</math> (так как матрица сингулярных чисел <math>\Sigma</math> - диагональная, то достаточно хранить вектор сингулярных чисел (элементы <math>\sigma_i</math>), а также две матрицы <math>U, V</math> левых и правых сингулярных векторов).

=== Свойства алгоритма ===
Метод Якоби является самым медленным из имеющихся алгоритмов вычисления сингулярных чисел и сингулярных векторов матрицы. Тем не менее, для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. Кроме того, в методе Якоби заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, за счет этого метод достигает высокой относительной точности.

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов для метода Якоби является линейным, вычислительная мощность алгоритма также линейна.

Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности (число итераций зависит от входных данных и порогового значения).

== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===

=== Локальность данных и вычислений ===

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

=== Масштабируемость алгоритма и его реализации ===
Данный раздел будет дополнен позже в связи с продлением срока сдачи (до 15 ноября).

=== Динамические характеристики и эффективность реализации алгоритма ===

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

=== Существующие реализации алгоритма ===
Несмотря на то, что метод Якоби является самым медленным из методов вычисления сингулярных чисел и векторов матрицы, он реализован в некоторых пакетах, например, LAPACK,scaLAPACK, Intel MKL<ref>Intel® Math Kernel Library,https://software.intel.com/ru-ru/node/521153#CAF8CA08-1C69-48F8-B0C3-B2B27AFA92FA</ref>.

== Литература ==

Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов

2016-10-15T00:22:55Z

Светлана: /* Информационный граф */

{{algorithm
| name = Метод Якоби вычисления сингулярных чисел и векторов
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math></math>
| pf_width = <math></math>
| input_data = <math>n^2</math>
| output_data = <math>2n^2 + n</math>
}}
Авторы статьи: [[Участник:Почернина Елена|Почернина Елена]] (группа 601), [[Участница:Светлана|Костюкова Светлана]] (группа 601)

== Свойства и структура алгоритма ==

=== Общее описание алгоритма ===

'''Метод Якоби''' позволяет найти сингулярное разложение плотной матрицы <math>G</math> размера <math>n</math>×<math>n</math>. Такое разложение можно записать в виде <math>G = U\Sigma V^*</math>, где <math>U, V</math> - унитарные (''в вещественном случае - ортогональные'') матрицы размера <math>n</math>×<math>n</math>, <math>\Sigma</math> - диагональная матрица размера <math>n</math>×<math>n</math> с вещественными положительными числами на главной диагонали. Диагональные элементы матрицы <math>\Sigma</math> называются сингулярными числами матрицы <math>G</math>, а столбцы матриц <math>U,V</math> - левыми и правыми сингулярными векторами.

''Неявный''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref> метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице <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>. Поэтому алгоритм называется ''методом односторонних вращений''.

Метод Якоби для вычисления сингулярных чисел и векторов является самым медленным из имеющихся алгоритмов поиска сингулярных чисел и векторов, но тем не менее, интерес к нему сохраняется. Для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. В частности, метод Якоби вычисляет сингулярные числа матрицы <math>G</math> с высокой точностью, если G может быть представлена в виде <math>G = DX</math>, где <math>D</math> – диагональная матрица, а <math>X</math> – хорошо обусловлена. В этом случае заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, в то время как другие алгоритмы включают в себя приведение матрицы к двухдиагональной форме, из-за чего и теряют все верные разряды во всех сингулярных числах, кроме старшего.

=== Математическое описание алгоритма ===
Исходные данные: матрица <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> - матрица правых сингулярных векторов.

Основной алгоритм заключается в следующем<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref>:

Для всех элементов матрицы <math>A^{(i)} = (G^TG)^{(i)}</math>, находящихся вне главной диагонали <math>(a_{12}, a_{13},\dots, a_{1n}, a_{23}, a_{24},\dots, a_{2n},\dots, a_{n-1n})</math> проверяем условие <math>|A_{jk}|<\varepsilon</math>:
Если <math>|A_{jk}^{(i)}| \le \varepsilon</math>, то переходим к следующему элементу,
Если <math>|A_{jk}^{(i)}| > \varepsilon</math>, то нужно выполнить вращение в плоскости <math>(jk)</math>, после чего перейти к следующей итерации алгоритма.
Если на очередной итерации не было совершено ни одного вращения, значит матрица <math>A</math> стала достаточна близка к диагональной, и можно перейти к вычислению искомых значений:

<math>\sigma_i = ||G(:,i)||_2</math>;

<math>U = [u_1,\dots,u_n]</math>, где <math>u_i = G(:,i)/\sigma_i</math>;

<math>V = J</math>, где <math>J</math> - накопленное произведение вращений Якоби.

Одностороннее вращение в плоскости <math>(jk)</math> осуществляет процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''':

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

Здесь <math>R(j,k,\theta)</math> - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 244-245)</ref>.

                                                                                
<math>
\begin{matrix}
& & & & j & & k & & & \\
\end{matrix}
</math>

<math>
J_i = 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>a_{jk}</math> не слишком мал"''. Подходящим критерием останова является условие
<math>\begin{align}
|a_{jk}| \ge \varepsilon\sqrt{a_{jj}a_{kk}}
\end{align}</math>.

=== Вычислительное ядро алгоритма ===

Вычислительное ядро составляют множественные вычисления элементов матрицы <math>a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk}</math> и <math>a_{kk} = (G^TG)_{kk}</math>, а также вычисление промежуточной матрицы <math>GR(j, k, \theta)</math> путем умножения исходной матрицы на матрицу вращения Якоби в процессе выполнения процедуры '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''.

После приведения матрицы <math>G^TG</math> к виду, достаточно близкому к диагональному, вычисляются сингулярные числа:

<math>\sigma_i = ||G(:,i)||_2</math>

Матрица левых сингулярных векторов <math>U = [u_1,..., u_n]</math>:

<math>u_i = G(:,i)/\sigma_i</math>

И (если требуется) матрица правых сингулярных векторов <math>V</math> путем накопления произведения вращений Якоби:

<math>
V = J,
J = JR(j, k, \theta)
</math>

=== Макроструктура алгоритма ===
Основную часть метода составляет процедура односторонних вращений матрицы (процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''). Рассмотрим подробнее структуру этой процедуры.

- Вычисление <math>a_{jk}, a_{jj}, a_{kk} (a_{ij} = (G^TG)_{ij})</math>. В связи с тем, что нет необходимости каждый раз вычислять произведение <math>G^TG</math>, для вычисления каждого элемента требуется <math>n</math> операций умножения и <math>(n-1)</math> операций сложения.

- Вычисление <math>\tau, c, s, t</math> не представляют значительной вычислительной сложности.

- Одностороннее вращение матрицы (умножение входной матрицы на матрицу вращения) требует порядка <math>O(n)</math> операций.

В итоге, процедура '''''One-Sided-Jacobi-Rotation''''' требует <math>O(n)</math> операций.

После выполнения алгоритма односторонних вращений необходимо вычислить искомые матрицы <math>S, U, V</math> (матрица сингулярных чисел и матрицы левых и правых векторов).

- Вычисление сингулярных чисел <math>\sigma_i</math> требует <math>n</math> операций умножения.

- Вычисление <math>U</math> требует <math>O(n^2)</math> операций. Аналогично для матрицы <math>V</math>.

=== Схема реализации последовательного алгоритма ===
Метод можно описать следующим образом<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-263)</ref>:

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> - накопленное произведение вращений Якоби.

Метод Якоби для нахождения сингулярных чисел и векторов матрицы <math>G</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

=== Последовательная сложность алгоритма ===
Для вычисления сингулярных чисел и векторов матрицы порядка <math>n</math> в последовательном варианте на каждой итерации алгоритма вызывается процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''', сложность которой равна <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>j</math> — меняется в диапазоне от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения;
* <math>k</math> — меняется в диапазоне от <math>j+1</math> до <math>n</math>, принимая все целочисленные значения.

'''Вторая''' группа вершин соответствует вычислению значений <math>c</math> и <math>s</math> для матрицы вращения

'''Третья''' группа вершин соответствует вычислению произведений матриц <math>GR(k,j, \theta)</math> и <math>JR(k,j, \theta)</math>

<center>
<div class="thumb">
<div class="thumbinner" style="width:{{#expr: 3 * (400 + 35) + 4 * (3 - 1) + 8}}px">
<gallery widths=400px heights=400px>
File:JacobiSVG.png|Рисунок 1. Граф одной итерации алгоритма без отображения входных и выходных данных.
</gallery>
</div>
</div>
</center>

=== Ресурс параллелизма алгоритма ===

=== Входные и выходные данные алгоритма ===
'''Входные данные''': плотная матрица <math>G</math> (элементы <math>g_{ij}</math>).
Дополнительные ограничения:
* <math>A= G^TG</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

'''Объём входных данных''': <math>n^2</math>

'''Выходные данные''': матрица <math>\Sigma = diag(\sigma_{i})</math>, где <math>\sigma_{i}</math> - сингулярные числа, матрица <math>U</math> левых сингулярных векторов и матрица <math>V</math> правых сингулярных векторов.

'''Объём выходных данных''': <math>2n^2 + n</math> (так как матрица сингулярных чисел <math>\Sigma</math> - диагональная, то достаточно хранить вектор сингулярных чисел (элементы <math>\sigma_i</math>), а также две матрицы <math>U, V</math> левых и правых сингулярных векторов).

=== Свойства алгоритма ===
Метод Якоби является самым медленным из имеющихся алгоритмов вычисления сингулярных чисел и сингулярных векторов матрицы. Тем не менее, для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. Кроме того, в методе Якоби заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, за счет этого метод достигает высокой относительной точности.

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов для метода Якоби является линейным, вычислительная мощность алгоритма также линейна.

Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности (число итераций зависит от входных данных и порогового значения).

== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===

=== Локальность данных и вычислений ===

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

=== Масштабируемость алгоритма и его реализации ===

=== Динамические характеристики и эффективность реализации алгоритма ===

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

=== Существующие реализации алгоритма ===
Метод Якоби нахождения сингулярных значений и векторов реализован в немногих пакетах, например, LAPACK, Intel MKL. Это не значит, что алгоритм не представляет никакого интереса для изучения. Для некоторых типов матриц метод способен вычислять сингулярные числа и сингулярные векторы намного точнее, чем другие методы. Но в связи с тем, что алгоритм считается медленным, он не включен во многие известные пакеты.

== Литература ==

Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов

2016-10-15T00:22:22Z

Светлана: /* Информационный граф */

{{algorithm
| name = Метод Якоби вычисления сингулярных чисел и векторов
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math></math>
| pf_width = <math></math>
| input_data = <math>n^2</math>
| output_data = <math>2n^2 + n</math>
}}
Авторы статьи: [[Участник:Почернина Елена|Почернина Елена]] (группа 601), [[Участница:Светлана|Костюкова Светлана]] (группа 601)

== Свойства и структура алгоритма ==

=== Общее описание алгоритма ===

'''Метод Якоби''' позволяет найти сингулярное разложение плотной матрицы <math>G</math> размера <math>n</math>×<math>n</math>. Такое разложение можно записать в виде <math>G = U\Sigma V^*</math>, где <math>U, V</math> - унитарные (''в вещественном случае - ортогональные'') матрицы размера <math>n</math>×<math>n</math>, <math>\Sigma</math> - диагональная матрица размера <math>n</math>×<math>n</math> с вещественными положительными числами на главной диагонали. Диагональные элементы матрицы <math>\Sigma</math> называются сингулярными числами матрицы <math>G</math>, а столбцы матриц <math>U,V</math> - левыми и правыми сингулярными векторами.

''Неявный''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref> метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице <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>. Поэтому алгоритм называется ''методом односторонних вращений''.

Метод Якоби для вычисления сингулярных чисел и векторов является самым медленным из имеющихся алгоритмов поиска сингулярных чисел и векторов, но тем не менее, интерес к нему сохраняется. Для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. В частности, метод Якоби вычисляет сингулярные числа матрицы <math>G</math> с высокой точностью, если G может быть представлена в виде <math>G = DX</math>, где <math>D</math> – диагональная матрица, а <math>X</math> – хорошо обусловлена. В этом случае заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, в то время как другие алгоритмы включают в себя приведение матрицы к двухдиагональной форме, из-за чего и теряют все верные разряды во всех сингулярных числах, кроме старшего.

=== Математическое описание алгоритма ===
Исходные данные: матрица <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> - матрица правых сингулярных векторов.

Основной алгоритм заключается в следующем<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref>:

Для всех элементов матрицы <math>A^{(i)} = (G^TG)^{(i)}</math>, находящихся вне главной диагонали <math>(a_{12}, a_{13},\dots, a_{1n}, a_{23}, a_{24},\dots, a_{2n},\dots, a_{n-1n})</math> проверяем условие <math>|A_{jk}|<\varepsilon</math>:
Если <math>|A_{jk}^{(i)}| \le \varepsilon</math>, то переходим к следующему элементу,
Если <math>|A_{jk}^{(i)}| > \varepsilon</math>, то нужно выполнить вращение в плоскости <math>(jk)</math>, после чего перейти к следующей итерации алгоритма.
Если на очередной итерации не было совершено ни одного вращения, значит матрица <math>A</math> стала достаточна близка к диагональной, и можно перейти к вычислению искомых значений:

<math>\sigma_i = ||G(:,i)||_2</math>;

<math>U = [u_1,\dots,u_n]</math>, где <math>u_i = G(:,i)/\sigma_i</math>;

<math>V = J</math>, где <math>J</math> - накопленное произведение вращений Якоби.

Одностороннее вращение в плоскости <math>(jk)</math> осуществляет процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''':

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

Здесь <math>R(j,k,\theta)</math> - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 244-245)</ref>.

                                                                                
<math>
\begin{matrix}
& & & & j & & k & & & \\
\end{matrix}
</math>

<math>
J_i = 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>a_{jk}</math> не слишком мал"''. Подходящим критерием останова является условие
<math>\begin{align}
|a_{jk}| \ge \varepsilon\sqrt{a_{jj}a_{kk}}
\end{align}</math>.

=== Вычислительное ядро алгоритма ===

Вычислительное ядро составляют множественные вычисления элементов матрицы <math>a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk}</math> и <math>a_{kk} = (G^TG)_{kk}</math>, а также вычисление промежуточной матрицы <math>GR(j, k, \theta)</math> путем умножения исходной матрицы на матрицу вращения Якоби в процессе выполнения процедуры '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''.

После приведения матрицы <math>G^TG</math> к виду, достаточно близкому к диагональному, вычисляются сингулярные числа:

<math>\sigma_i = ||G(:,i)||_2</math>

Матрица левых сингулярных векторов <math>U = [u_1,..., u_n]</math>:

<math>u_i = G(:,i)/\sigma_i</math>

И (если требуется) матрица правых сингулярных векторов <math>V</math> путем накопления произведения вращений Якоби:

<math>
V = J,
J = JR(j, k, \theta)
</math>

=== Макроструктура алгоритма ===
Основную часть метода составляет процедура односторонних вращений матрицы (процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''). Рассмотрим подробнее структуру этой процедуры.

- Вычисление <math>a_{jk}, a_{jj}, a_{kk} (a_{ij} = (G^TG)_{ij})</math>. В связи с тем, что нет необходимости каждый раз вычислять произведение <math>G^TG</math>, для вычисления каждого элемента требуется <math>n</math> операций умножения и <math>(n-1)</math> операций сложения.

- Вычисление <math>\tau, c, s, t</math> не представляют значительной вычислительной сложности.

- Одностороннее вращение матрицы (умножение входной матрицы на матрицу вращения) требует порядка <math>O(n)</math> операций.

В итоге, процедура '''''One-Sided-Jacobi-Rotation''''' требует <math>O(n)</math> операций.

После выполнения алгоритма односторонних вращений необходимо вычислить искомые матрицы <math>S, U, V</math> (матрица сингулярных чисел и матрицы левых и правых векторов).

- Вычисление сингулярных чисел <math>\sigma_i</math> требует <math>n</math> операций умножения.

- Вычисление <math>U</math> требует <math>O(n^2)</math> операций. Аналогично для матрицы <math>V</math>.

=== Схема реализации последовательного алгоритма ===
Метод можно описать следующим образом<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-263)</ref>:

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> - накопленное произведение вращений Якоби.

Метод Якоби для нахождения сингулярных чисел и векторов матрицы <math>G</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

=== Последовательная сложность алгоритма ===
Для вычисления сингулярных чисел и векторов матрицы порядка <math>n</math> в последовательном варианте на каждой итерации алгоритма вызывается процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''', сложность которой равна <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>j</math> — меняется в диапазоне от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения;
* <math>k</math> — меняется в диапазоне от <math>j+1</math> до <math>n</math>, принимая все целочисленные значения.

'''Вторая''' группа вершин соответствует вычислению значений <math>c</math> и <math>s</math> для матрицы вращения

'''Третья''' группа вершин соответствует вычислению произведений матриц <math>GR(k,j, \theta)</math> и <math>JR(k,j, \theta)</math>

<center>
<div class="thumb">
<div class="thumbinner" style="width:{{#expr: 3 * (400 + 35) + 4 * (3 - 1) + 8}}px">
<gallery widths=400px heights=400px>
File:JacobiSVG.png|Рисунок 1. Граф одной итерации алгоритма без отображения входных и выходных данных.</math>.
</gallery>
</div>
</div>
</center>

=== Ресурс параллелизма алгоритма ===

=== Входные и выходные данные алгоритма ===
'''Входные данные''': плотная матрица <math>G</math> (элементы <math>g_{ij}</math>).
Дополнительные ограничения:
* <math>A= G^TG</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

'''Объём входных данных''': <math>n^2</math>

'''Выходные данные''': матрица <math>\Sigma = diag(\sigma_{i})</math>, где <math>\sigma_{i}</math> - сингулярные числа, матрица <math>U</math> левых сингулярных векторов и матрица <math>V</math> правых сингулярных векторов.

'''Объём выходных данных''': <math>2n^2 + n</math> (так как матрица сингулярных чисел <math>\Sigma</math> - диагональная, то достаточно хранить вектор сингулярных чисел (элементы <math>\sigma_i</math>), а также две матрицы <math>U, V</math> левых и правых сингулярных векторов).

=== Свойства алгоритма ===
Метод Якоби является самым медленным из имеющихся алгоритмов вычисления сингулярных чисел и сингулярных векторов матрицы. Тем не менее, для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. Кроме того, в методе Якоби заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, за счет этого метод достигает высокой относительной точности.

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов для метода Якоби является линейным, вычислительная мощность алгоритма также линейна.

Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности (число итераций зависит от входных данных и порогового значения).

== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===

=== Локальность данных и вычислений ===

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

=== Масштабируемость алгоритма и его реализации ===

=== Динамические характеристики и эффективность реализации алгоритма ===

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

=== Существующие реализации алгоритма ===
Метод Якоби нахождения сингулярных значений и векторов реализован в немногих пакетах, например, LAPACK, Intel MKL. Это не значит, что алгоритм не представляет никакого интереса для изучения. Для некоторых типов матриц метод способен вычислять сингулярные числа и сингулярные векторы намного точнее, чем другие методы. Но в связи с тем, что алгоритм считается медленным, он не включен во многие известные пакеты.

== Литература ==

Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов

2016-10-14T14:22:23Z

Светлана: /* Вычислительное ядро алгоритма */

{{algorithm
| name = Метод Якоби вычисления сингулярных чисел и векторов
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math></math>
| pf_width = <math></math>
| input_data = <math>n^2</math>
| output_data = <math>2n^2 + n</math>
}}
Авторы статьи: [[Участник:Почернина Елена|Почернина Елена]] (группа 601), [[Участница:Светлана|Костюкова Светлана]] (группа 601)

== Свойства и структура алгоритма ==

=== Общее описание алгоритма ===

'''Метод Якоби''' позволяет найти сингулярное разложение плотной матрицы <math>G</math> размера <math>n</math>×<math>n</math>. Такое разложение можно записать в виде <math>G = U\Sigma V^*</math>, где <math>U, V</math> - унитарные (''в вещественном случае - ортогональные'') матрицы размера <math>n</math>×<math>n</math>, <math>\Sigma</math> - диагональная матрица размера <math>n</math>×<math>n</math> с вещественными положительными числами на главной диагонали. Диагональные элементы матрицы <math>\Sigma</math> называются сингулярными числами матрицы <math>G</math>, а столбцы матриц <math>U,V</math> - левыми и правыми сингулярными векторами.

''Неявный'' метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице <math>A = G^TG</math>.
Это значит, что на каждом шаге вычисляется вращение Якоби <math>J</math>, с помощью которого матрица <math>G^TG</math> неявно пересчитывется в <math>J^TG^TGJ</math>

Метод Якоби для вычисления сингулярных чисел и векторов является самым медленным из имеющихся алгоритмов поиска сингулярных чисел и векторов, но тем не менее, интерес к нему сохраняется. Для некоторых некоторых типов матриц <math>G</math> он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. В частности, метод Якоби вычисляет сингулярные числа матрицы <math>G</math> с высокой точностью, если G может быть представлена в виде <math>G = DX</math>, где <math>D</math> – диагональная матрица, а <math>X</math> – хорошо обусловлена. В этом случае заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, в то время как другие алгоритмы включают в себя приведение матрицы к двухдиагональной форме, из-за чего и теряют все верные разряды во всех сингулярных числах, кроме старшего.

=== Математическое описание алгоритма ===
Исходные данные: матрица <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> - матрица правых сингулярных векторов.

Основной алгоритм заключается в следующем<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)</ref>:

Для всех элементов матрицы <math>A^{(i)} = (G^TG)^{(i)}</math>, находящихся вне главной диагонали <math>(a_{12}, a_{13},\dots, a_{1n}, a_{23}, a_{24},\dots, a_{2n},\dots, a_{n-1n})</math> проверяем условие <math>|A_{jk}|<\varepsilon</math>:
Если <math>|A_{jk}^{(i)}| \le \varepsilon</math>, то переходим к следующему элементу,
Если <math>|A_{jk}^{(i)}| > \varepsilon</math>, то нужно выполнить вращение в плоскости <math>(jk)</math>, после чего перейти к следующей итерации алгоритма.
Если на очередной итерации не было совершено ни одного вращения, значит матрица <math>A</math> стала достаточна близка к диагональной, и можно перейти к вычислению искомых значений:

<math>\sigma_i = ||G(:,i)||_2</math>;

<math>U = [u_1,\dots,u_n]</math>, где <math>u_i = G(:,i)/\sigma_i</math>;

<math>V = J</math>, где <math>J</math> - накопленное произведение вращений Якоби.

Одностороннее вращение в плоскости <math>(jk)</math> осуществляет процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>''''':

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

Здесь <math>R(j,k,\theta)</math> - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 244-245)</ref>.

                                                                                
<math>
\begin{matrix}
& & & & j & & k & & & \\
\end{matrix}
</math>

<math>
J_i = 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>a_{jk}</math> не слишком мал"''. Подходящим критерием останова является условие
<math>\begin{align}
|a_{jk}| \ge \varepsilon\sqrt{a_{jj}a_{kk}}
\end{align}</math>.

=== Вычислительное ядро алгоритма ===

Вычислительное ядро составляют множественные вычисления элементов матрицы <math>a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk}</math> и <math>a_{kk} = (G^TG)_{kk}</math>, а также вычисление промежуточной матрицы <math>GR(j, k, \theta)</math> путем умножения исходной матрицы на матрицу вращения Якоби в процессе выполнения процедуры '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''.

После приведения матрицы <math>G^TG</math> к виду, достаточно близкому к диагональному, вычисляются сингулярные числа:

<math>\sigma_i = ||G(:,i)||_2</math>

Матрица левых сингулярных векторов <math>U = [u_1,..., u_n]</math>:

<math>u_i = G(:,i)/\sigma_i</math>

И (если требуется) матрица правых сингулярных векторов <math>V</math> путем накопления произведения вращений Якоби:

<math>
V = J,
J = JR(j, k, \theta)
</math>

=== Макроструктура алгоритма ===
Основную часть метода составляет процедура односторонних вращений матрицы (процедура '''''One-Sided-Jacobi-Rotation<math>(G,j,k)</math>'''''). Рассмотрим подробнее структуру этой процедуры.

- Вычисление <math>a_{jk}, a_{jj}, a_{kk} (a_{ij} = (G^TG)_{ij})</math>. В связи с тем, что нет необходимости каждый раз вычислять произведение <math>G^TG</math>, для вычисления каждого элемента требуется <math>n</math> операций умножения и <math>(n-1)</math> операций сложения.

- Вычисление <math>\tau, c, s, t</math> не представляют значительной вычислительной сложности.

- Одностороннее вращение матрицы (умножение входной матрицы на матрицу вращения) требует порядка <math>O(n)</math> операций.

В итоге, процедура '''''One-Sided-Jacobi-Rotation''''' требует <math>O(n)</math> операций.

После выполнения алгоритма односторонних вращений необходимо вычислить искомые матрицы <math>S, U, V</math> (матрица сингулярных чисел и матрицы левых и правых векторов).

- Вычисление сингулярных чисел <math>\sigma_i</math> требует <math>n</math> операций умножения.

- Вычисление <math>U</math> требует <math>O(n^2)</math> операций. Аналогично для матрицы <math>V</math>.

=== Схема реализации последовательного алгоритма ===
Метод можно описать следующим образом<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-263)</ref>:

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> - накопленное произведение вращений Якоби.

Метод Якоби для нахождения сингулярных чисел и векторов матрицы <math>G</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

=== Последовательная сложность алгоритма ===
Для вычисления сингулярных чисел и векторов матрицы порядка n в последовательном варианте требуется:

* <math></math> делений,
* <math></math> сложений (вычитаний),
* <math></math> умножений.

Метод Якоби вычисления сингулярных чисел и векторов при классификации по последовательной сложности относится к алгоритмам с ''кубической сложностью''.

=== Информационный граф ===
Опишем граф алгоритма как аналитически, так и в виде рисунка.

'''Первая''' группа вершин расположена в двумерной области, соответствующая ей операция - вычисление <math>a_{jj}, a_{jk}, a_{kk}</math>.
Естественно введённые координаты области таковы:
* <math>j</math> — меняется в диапазоне от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения;
* <math>k</math> — меняется в диапазоне от <math>j+1</math> до <math>n</math>, принимая все целочисленные значения.

=== Ресурс параллелизма алгоритма ===

=== Входные и выходные данные алгоритма ===
'''Входные данные''': плотная матрица <math>G</math> (элементы <math>g_{ij}</math>).
Дополнительные ограничения:
* <math>A= G^TG</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

'''Объём входных данных''': <math>n^2</math>

'''Выходные данные''': матрица <math>\Sigma = diag(\sigma_{i})</math>, где <math>\sigma_{i}</math> - сингулярные числа, матрица <math>U</math> левых сингулярных векторов и матрица <math>V</math> правых сингулярных векторов.

'''Объём выходных данных''': <math>2n^2 + n</math> (так как матрица сингулярных чисел <math>\Sigma</math> - диагональная, то достаточно хранить вектор сингулярных чисел (элементы <math>\sigma_i</math>), а также две матрицы <math>U, V</math> левых и правых сингулярных векторов).

=== Свойства алгоритма ===

== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===

=== Локальность данных и вычислений ===

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

=== Масштабируемость алгоритма и его реализации ===

=== Динамические характеристики и эффективность реализации алгоритма ===

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

=== Существующие реализации алгоритма ===
Метод Якоби нахождения сингулярных значений и векторов реализован в немногих пакетах, например, LAPACK, Intel MKL. Это не значит, что алгоритм не представляет никакого интереса для изучения. Для некоторых типов матриц метод способен вычислять сингулярные числа и сингулярные векторы намного точнее, чем другие методы. Но в связи с тем, что алгоритм считается медленным, он не включен во многие известные пакеты.

== Литература ==

Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов

2016-10-13T15:30:39Z

Светлана: /* Вычислительное ядро алгоритма */

Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов

2016-10-13T14:21:10Z

Светлана: /* Математическое описание алгоритма */

Обсуждение участника:Lineprinter

2016-10-11T22:46:11Z

Светлана: /* Перенос */

== Перенос ==

Статья: algowiki-project.org/ru/Решение_задачи_Коши_для_системы_ОДУ_методом_Рунге-Кутты_4_порядка;
Выполняют: Алимов Дамир и Антон Лесничий

Статья: Алгоритм кластеризации, основанный на минимальном остовном дереве;
Выполняет: Стандрик Антон Сергеевич (Astandrik)

Статья: [[Нахождение_собственных_чисел_квадратной_матрицы_методом_QR_разложения_(3)|Нахождение собственных чисел квадратной матрицы методом QR разложения (3)]];
Выполняют: [[Участник:Смирнова_Александра|Смирнова А.С.]], [[Участник:Александра_Киямова|Киямова А.]]

Статья: [[Метод_Якоби_вычисления_сингулярных_чисел_и_векторов]];
Выполняют: [[Участница:Почернина_Елена| Почернина Елена]] и [[Участница:Светлана|Костюкова Светлана]]

Обсуждение участника:Lineprinter

2016-10-11T22:33:52Z

Светлана: /* Перенос */

== Перенос ==

Статья: algowiki-project.org/ru/Решение_задачи_Коши_для_системы_ОДУ_методом_Рунге-Кутты_4_порядка;
Выполняют: Алимов Дамир и Антон Лесничий

Статья: Алгоритм кластеризации, основанный на минимальном остовном дереве;
Выполняет: Стандрик Антон Сергеевич (Astandrik)

Статья: [[Нахождение_собственных_чисел_квадратной_матрицы_методом_QR_разложения_(3)|Нахождение собственных чисел квадратной матрицы методом QR разложения (3)]];
Выполняют: [[Участник:Смирнова_Александра|Смирнова А.С.]], [[Участник:Александра_Киямова|Киямова А.]]

Статья: [[Метод_Якоби_вычисления_сингулярных_чисел_и_векторов]];
Выполняют: Почернина Елена и Костюкова Светлана