Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов: различия между версиями
Строка 24: | Строка 24: | ||
Основной алгоритм заключается в следующем: | Основной алгоритм заключается в следующем: | ||
+ | |||
Для всех элементов матрицы <math>A^{(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}|<\epsilon</math>: | Для всех элементов матрицы <math>A^{(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}|<\epsilon</math>: | ||
Если <math>|A_{jk}^{(i)}| \le \epsilon</math>, то переходим к следующему элементу, | Если <math>|A_{jk}^{(i)}| \le \epsilon</math>, то переходим к следующему элементу, | ||
Строка 34: | Строка 35: | ||
<math>V = J</math>, где <math>J</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> - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов. | ||
+ | |||
+ | | ||
+ | <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> | ||
=== Вычислительное ядро алгоритма === | === Вычислительное ядро алгоритма === |
Версия 12:58, 14 октября 2016
Метод Якоби вычисления сингулярных чисел и векторов | |
Последовательный алгоритм | |
Последовательная сложность | [math]O(n^3)[/math] |
Объём входных данных | [math][/math] |
Объём выходных данных | [math][/math] |
Параллельный алгоритм | |
Высота ярусно-параллельной формы | [math][/math] |
Ширина ярусно-параллельной формы | [math][/math] |
Авторы статьи: Почернина Елена (группа 601), Костюкова Светлана (группа 601)
Содержание
- 1 Свойства и структура алгоритма
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 3 Литература
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
Метод Якоби для вычисления сингулярных чисел и векторов является самым медленным из имеющихся алгоритмов поиска сингулярных чисел и векторов, но тем не менее, интерес к нему сохраняется. Для некоторых некоторых типов матриц [math]G[/math] он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. В частности, метод Якоби вычисляет сингулярные числа матрицы [math]G[/math] с высокой точностью, если G может быть представлена в виде [math]G = DX[/math], где [math]D[/math] – диагональная матрица, а [math]X[/math] – хорошо обусловлена. В этом случае заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, в то время как другие алгоритмы включают в себя приведение матрицы к двухдиагональной форме, из-за чего и теряют все верные разряды во всех сингулярных числах, кроме старшего.
1.2 Математическое описание алгоритма
Исходные данные: матрица [math]G[/math] (элементы [math]g_{ij}, i, j = 1, \ldots, n[/math]).
Вычисляемые данные: матрица [math]\Sigma = diag(\sigma_{i})[/math], где [math]\sigma_{i}[/math] - сингулярные числа, [math]U[/math] - матрица левых сингулярных векторов, [math]V[/math] - матрица правых сингулярных векторов.
Матрица [math]A = G^TG[/math].
Основной алгоритм заключается в следующем:
Для всех элементов матрицы [math]A^{(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}|\lt \epsilon[/math]: Если [math]|A_{jk}^{(i)}| \le \epsilon[/math], то переходим к следующему элементу, Если [math]|A_{jk}^{(i)}| \gt \epsilon[/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] - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов.
[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]
1.3 Вычислительное ядро алгоритма
Вычислительное ядро составляют множественные вычисления элементов матрицы [math]a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk}[/math] и [math]a_{kk} = (G^TG)_{kk}[/math], а также вычисление промежуточной матрицы [math]GJ[/math].
1.4 Макроструктура алгоритма
1.5 Схема реализации последовательного алгоритма
Метод можно описать следующим образом:
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
1.6 Последовательная сложность алгоритма
Для вычисления сингулярных чисел и векторов матрицы порядка n в последовательном варианте требуется:
- [math][/math] делений,
- [math][/math] сложений (вычитаний),
- [math][/math] умножений.
Метод Якоби вычисления сингулярных чисел и векторов при классификации по последовательной сложности относится к алгоритмам с кубической сложностью.
1.7 Информационный граф
Опишем граф алгоритма как аналитически, так и в виде рисунка.
Первая группа вершин расположена в двумерной области, соответствующая ей операция - вычисление [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], принимая все целочисленные значения.
1.8 Ресурс параллелизма алгоритма
1.9 Входные и выходные данные алгоритма
Входные данные: плотная матрица [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]
1.10 Свойства алгоритма
2 Программная реализация алгоритма
2.1 Масштабируемость алгоритма и его реализации
2.2 Существующие реализации алгоритма
Метод Якоби нахождения сингулярных значений и векторов реализован в пакете LAPACK. В связи с тем, что алгоритм считается медленным, он не включен во многие известные пакеты.
3 Литература
1. Дж. Деммель Вычислительная линейная алгебра. Изд. Мир, 2001.