Участник:Почернина Елена/Метод Якоби вычисления сингулярных чисел и векторов
| Метод Якоби вычисления сингулярных чисел и векторов | |
| Последовательный алгоритм | |
| Последовательная сложность | [math]\displaystyle{ O(n^3) }[/math] |
| Объём входных данных | [math]\displaystyle{ n^2 }[/math] |
| Объём выходных данных | [math]\displaystyle{ 2n^2 + n }[/math] |
| Параллельный алгоритм | |
| Высота ярусно-параллельной формы | [math]\displaystyle{ O(n^2) }[/math] |
| Ширина ярусно-параллельной формы | [math]\displaystyle{ O(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]\displaystyle{ G }[/math] размера [math]\displaystyle{ n }[/math]×[math]\displaystyle{ n }[/math]. Такое разложение можно записать в виде [math]\displaystyle{ G = U\Sigma V^* }[/math], где [math]\displaystyle{ U, V }[/math] - унитарные (в вещественном случае - ортогональные) матрицы размера [math]\displaystyle{ n }[/math]×[math]\displaystyle{ n }[/math], [math]\displaystyle{ \Sigma }[/math] - диагональная матрица размера [math]\displaystyle{ n }[/math]×[math]\displaystyle{ n }[/math] с вещественными положительными числами на главной диагонали. Диагональные элементы матрицы [math]\displaystyle{ \Sigma }[/math] называются сингулярными числами матрицы [math]\displaystyle{ G }[/math], а столбцы матриц [math]\displaystyle{ U,V }[/math] - левыми и правыми сингулярными векторами.
Неявный[1] метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице [math]\displaystyle{ A = G^TG }[/math]. Это значит, что на каждом шаге вычисляется вращение Якоби [math]\displaystyle{ J }[/math], с помощью которого матрица [math]\displaystyle{ G^TG }[/math] неявно пересчитывется в [math]\displaystyle{ J^TG^TGJ }[/math]; вращение выбрано так, чтобы пара внедиагональных элементов из [math]\displaystyle{ G^TG }[/math] обратилась в нули в матрице [math]\displaystyle{ J^TG^TGJ }[/math]. При этом ни [math]\displaystyle{ G^TG }[/math], ни [math]\displaystyle{ J^TG^TGJ }[/math] не вычисляются в явном виде; вместо них вычисляется матрица [math]\displaystyle{ GJ }[/math]. Поэтому алгоритм называется методом односторонних вращений.
Метод Якоби для вычисления сингулярных чисел и векторов является самым медленным из имеющихся алгоритмов поиска сингулярных чисел и векторов, но тем не менее, интерес к нему сохраняется. Для некоторых некоторых типов матриц [math]\displaystyle{ G }[/math] он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. В частности, метод Якоби вычисляет сингулярные числа матрицы [math]\displaystyle{ G }[/math] с высокой точностью, если G может быть представлена в виде [math]\displaystyle{ G = DX }[/math], где [math]\displaystyle{ D }[/math] – диагональная матрица, а [math]\displaystyle{ X }[/math] – хорошо обусловлена. В этом случае заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, в то время как другие алгоритмы включают в себя приведение матрицы к двухдиагональной форме, из-за чего и теряют все верные разряды во всех сингулярных числах, кроме старшего.
Математическое описание алгоритма
Исходные данные: матрица [math]\displaystyle{ G }[/math] (элементы [math]\displaystyle{ g_{ij}, i, j = 1, \ldots, n }[/math]).
Вычисляемые данные: матрица [math]\displaystyle{ \Sigma = diag(\sigma_{i}) }[/math], где [math]\displaystyle{ \sigma_{i} }[/math] - сингулярные числа, [math]\displaystyle{ U }[/math] - матрица левых сингулярных векторов, [math]\displaystyle{ V }[/math] - матрица правых сингулярных векторов.
Основной алгоритм заключается в следующем[2]:
Для всех элементов матрицы [math]\displaystyle{ A^{(i)} = (G^TG)^{(i)} }[/math], находящихся вне главной диагонали [math]\displaystyle{ (a_{12}, a_{13},\dots, a_{1n}, a_{23}, a_{24},\dots, a_{2n},\dots, a_{n-1n}) }[/math] проверяем условие [math]\displaystyle{ |A_{jk}|\lt \varepsilon }[/math]: Если [math]\displaystyle{ |A_{jk}^{(i)}| \le \varepsilon }[/math], то переходим к следующему элементу, Если [math]\displaystyle{ |A_{jk}^{(i)}| \gt \varepsilon }[/math], то нужно выполнить вращение в плоскости [math]\displaystyle{ (j,k) }[/math], после чего перейти к следующей итерации алгоритма. Если на очередной итерации не было совершено ни одного вращения, значит матрица [math]\displaystyle{ A }[/math] стала достаточна близка к диагональной, и можно перейти к вычислению искомых значений:
[math]\displaystyle{ \sigma_i = ||G(:,i)||_2 }[/math];
[math]\displaystyle{ U = [u_1,\dots,u_n] }[/math], где [math]\displaystyle{ u_i = G(:,i)/\sigma_i }[/math];
[math]\displaystyle{ V = J }[/math], где [math]\displaystyle{ J }[/math] - накопленное произведение вращений Якоби.
Одностороннее вращение в плоскости [math]\displaystyle{ (j,k) }[/math] осуществляет процедура One-Sided-Jacobi-Rotation[math]\displaystyle{ (G,j,k) }[/math]:
proc One-Sided-Jacobi-Rotation[math]\displaystyle{ (G,j,k) }[/math] вычислить [math]\displaystyle{ a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk} }[/math] и [math]\displaystyle{ a_{kk} = (G^TG)_{kk} }[/math] if [math]\displaystyle{ |a_{jk}| }[/math] не слишком мал [math]\displaystyle{ \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]\displaystyle{ J = JR(j,k,\theta) }[/math] end if end if
Здесь [math]\displaystyle{ R(j,k,\theta) }[/math] - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов[3].
[math]\displaystyle{ \begin{matrix} & & & & j & & k & & & \\ \end{matrix} }[/math]
[math]\displaystyle{ J_i = R(j,k,\theta) = \begin{matrix} \\ \\ \\ j \\ \\ k \\ \\ \\ \\ \end{matrix} }[/math] [math]\displaystyle{ \begin{bmatrix} & 1 & & & & & & & & \\ & & 1 & & & & & & & \\ & & & \ddots & & & & & & \\ & & & & \cos(\theta) & & -\sin(\theta) & & & \\ & & & & & \ddots & & & & \\ & & & & \sin(\theta) & & \cos(\theta) & & & \\ & & & & & & & \ddots & & \\ & & & & & & & & 1 & \\ & & & & & & & & & 1 \\ \end{bmatrix} }[/math]
Для завершения описания алгоритма, нужно указать критерий останова, то есть объяснить смысл условия "если [math]\displaystyle{ a_{jk} }[/math] не слишком мал". Подходящим критерием останова является условие
[math]\displaystyle{ \begin{align}
|a_{jk}| \ge \varepsilon\sqrt{a_{jj}a_{kk}}
\end{align} }[/math].
Вычислительное ядро алгоритма
Вычислительное ядро составляют множественные вычисления элементов матрицы [math]\displaystyle{ a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk} }[/math] и [math]\displaystyle{ a_{kk} = (G^TG)_{kk} }[/math], а также вычисление промежуточной матрицы [math]\displaystyle{ GR(j, k, \theta) }[/math] путем умножения исходной матрицы на матрицу вращения Якоби в процессе выполнения процедуры One-Sided-Jacobi-Rotation[math]\displaystyle{ (G,j,k) }[/math].
После приведения матрицы [math]\displaystyle{ G^TG }[/math] к виду, достаточно близкому к диагональному, вычисляются сингулярные числа:
[math]\displaystyle{ \sigma_i = ||G(:,i)||_2 }[/math]
Матрица левых сингулярных векторов [math]\displaystyle{ U = [u_1,..., u_n] }[/math]:
[math]\displaystyle{ u_i = G(:,i)/\sigma_i }[/math]
И (если требуется) матрица правых сингулярных векторов [math]\displaystyle{ V }[/math] путем накопления произведения вращений Якоби:
[math]\displaystyle{ V = J, J = JR(j, k, \theta) }[/math]
Макроструктура алгоритма
Основную часть метода составляет процедура односторонних вращений матрицы (процедура One-Sided-Jacobi-Rotation[math]\displaystyle{ (G,j,k) }[/math]). Рассмотрим подробнее структуру этой процедуры.
- Вычисление [math]\displaystyle{ a_{jk}, a_{jj}, a_{kk} (a_{ij} = (G^TG)_{ij}) }[/math]. В связи с тем, что нет необходимости каждый раз вычислять произведение [math]\displaystyle{ G^TG }[/math], для вычисления каждого элемента требуется [math]\displaystyle{ n }[/math] операций умножения и [math]\displaystyle{ (n-1) }[/math] операций сложения.
- Вычисление [math]\displaystyle{ \tau, c, s, t }[/math] не представляют значительной вычислительной сложности.
- Одностороннее вращение матрицы (умножение входной матрицы на матрицу вращения) требует порядка [math]\displaystyle{ O(n) }[/math] операций.
В итоге, процедура One-Sided-Jacobi-Rotation требует [math]\displaystyle{ O(n) }[/math] операций.
После выполнения алгоритма односторонних вращений необходимо вычислить искомые матрицы [math]\displaystyle{ S, U, V }[/math] (матрица сингулярных чисел и матрицы левых и правых векторов).
- Вычисление сингулярных чисел [math]\displaystyle{ \sigma_i }[/math] требует [math]\displaystyle{ n }[/math] операций умножения.
- Вычисление [math]\displaystyle{ U }[/math] требует [math]\displaystyle{ O(n^2) }[/math] операций. Аналогично для матрицы [math]\displaystyle{ V }[/math].
Схема реализации последовательного алгоритма
Метод можно описать следующим образом[4]:
repeat
for [math]\displaystyle{ j=1 }[/math] to [math]\displaystyle{ n-1 }[/math]
for [math]\displaystyle{ k=j+1 }[/math] to [math]\displaystyle{ n }[/math]
обратиться к процедуре One-Sided-Jacobi-Rotation[math]\displaystyle{ (G,j,k) }[/math]
end for
end for
пока [math]\displaystyle{ G^TG }[/math] не станет достаточно близка к диагональной матрице
Положить [math]\displaystyle{ \sigma_{i} = ||G(:,i)||_2 }[/math] (норма [math]\displaystyle{ i }[/math]-го столбца в [math]\displaystyle{ G }[/math])
Положить [math]\displaystyle{ U = [u_{1},\dots, u_{n}] }[/math], где [math]\displaystyle{ u_{i} = G(:,i)/\sigma_{i} }[/math]
Положить [math]\displaystyle{ V = J }[/math], где [math]\displaystyle{ J }[/math] - накопленное произведение вращений Якоби.
Метод Якоби для нахождения сингулярных чисел и векторов матрицы [math]\displaystyle{ G }[/math] также использует метод односторонних вращений (процедура One-Sided-Jacobi-Rotation). Алгоритм выглядит следующим образом:
proc One-Sided-Jacobi-Rotation[math]\displaystyle{ (G,j,k) }[/math] вычислить [math]\displaystyle{ a_{jj} = (G^TG)_{jj}, a_{jk} = (G^TG)_{jk} }[/math] и [math]\displaystyle{ a_{kk} = (G^TG)_{kk} }[/math] if [math]\displaystyle{ |a_{jk}| }[/math] не слишком мал [math]\displaystyle{ \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]\displaystyle{ J = JR(j,k,\theta) }[/math] end if end if
Последовательная сложность алгоритма
Для вычисления сингулярных чисел и векторов матрицы порядка [math]\displaystyle{ n }[/math] в последовательном варианте на каждой итерации алгоритма вызывается процедура One-Sided-Jacobi-Rotation[math]\displaystyle{ (G,j,k) }[/math], сложность которой равна [math]\displaystyle{ O(n) }[/math]. Так как выбор индексов [math]\displaystyle{ j }[/math] и [math]\displaystyle{ k }[/math] осуществляется путем перебора внедиагональных элементов матрицы, для полного прохода потребуется [math]\displaystyle{ \frac{n(n-1)}{2} }[/math] итераций.
Таким образом, метод Якоби вычисления сингулярных чисел и векторов при классификации по последовательной сложности относится к алгоритмам с кубической сложностью.
Информационный граф
Граф алгоритма состоит из трех групп вершин:
Первая группа вершин расположена в двумерной области, соответствующая ей операция - вычисление [math]\displaystyle{ a_{jj}, a_{jk}, a_{kk} }[/math]. Естественно введённые координаты области таковы:
- [math]\displaystyle{ j }[/math] — меняется в диапазоне от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ n-1 }[/math], принимая все целочисленные значения;
- [math]\displaystyle{ k }[/math] — меняется в диапазоне от [math]\displaystyle{ j+1 }[/math] до [math]\displaystyle{ n }[/math], принимая все целочисленные значения.
Вторая группа вершин соответствует вычислению значений [math]\displaystyle{ c }[/math] и [math]\displaystyle{ s }[/math] для матрицы вращения
Третья группа вершин соответствует вычислению произведений матриц [math]\displaystyle{ GR(k,j, \theta) }[/math] и [math]\displaystyle{ JR(k,j, \theta) }[/math]
-
Рисунок 1. Граф одной итерации алгоритма без отображения входных и выходных данных. -
![Рисунок 2. Граф вычисления значений элементов [math]\displaystyle{ a_{jj}, a_{jk}, a_{kk} }[/math].](/w/pool/images/thumb/a/a4/Akj.png/400px-Akj.png)
Рисунок 2. Граф вычисления значений элементов [math]\displaystyle{ a_{jj}, a_{jk}, a_{kk} }[/math].
![Рисунок 2. Граф вычисления значений элементов [math]\displaystyle{ a_{jj}, a_{jk}, a_{kk} }[/math].](/ru/%D0%A4%D0%B0%D0%B9%D0%BB:Akj.png)
-
![Рисунок 3. Граф одной итерации алгоритма с вычислением матрицы вращения [math]\displaystyle{ R }[/math] отображением входных и выходных данных.](/w/pool/images/thumb/9/94/JacobiSVD2.png/400px-JacobiSVD2.png)
Рисунок 3. Граф одной итерации алгоритма с вычислением матрицы вращения [math]\displaystyle{ R }[/math] отображением входных и выходных данных.
![Рисунок 3. Граф одной итерации алгоритма с вычислением матрицы вращения [math]\displaystyle{ R }[/math] отображением входных и выходных данных.](/ru/%D0%A4%D0%B0%D0%B9%D0%BB:JacobiSVD2.png)
-
![Рисунок 4. Граф вычисления значений элементов матрицы [math]\displaystyle{ GR }[/math]](/w/pool/images/thumb/1/13/Gij.png/400px-Gij.png)
Рисунок 4. Граф вычисления значений элементов матрицы [math]\displaystyle{ GR }[/math] -
![Рисунок 5. Граф вычисления значений элементов матрицы [math]\displaystyle{ JR }[/math]](/w/pool/images/thumb/2/24/Jij.png/400px-Jij.png)
Рисунок 5. Граф вычисления значений элементов матрицы [math]\displaystyle{ JR }[/math]
![Рисунок 4. Граф вычисления значений элементов матрицы [math]\displaystyle{ GR }[/math]](/ru/%D0%A4%D0%B0%D0%B9%D0%BB:Gij.png)
![Рисунок 5. Граф вычисления значений элементов матрицы [math]\displaystyle{ JR }[/math]](/ru/%D0%A4%D0%B0%D0%B9%D0%BB:Jij.png)
Ресурс параллелизма алгоритма
В данном алгоритме существует возможность параллельного вычисления вращений Якоби[5][6]. Максимальное количество вращений в данном методе равно [math]\displaystyle{ \frac{n(n-1)}{2} }[/math]. В параллельной реализации, если [math]\displaystyle{ n }[/math] - четное, можно выполнить [math]\displaystyle{ \frac{n}{2} }[/math] вращений параллельно, если [math]\displaystyle{ n }[/math] - нечетное, то [math]\displaystyle{ \frac{n-1}{2} }[/math] вращений.
Основную сложность в процедуре вращения One-Sided-Jacobi-Rotation[math]\displaystyle{ (G,j,k) }[/math] представляет вычисление элементов [math]\displaystyle{ a_{jj}, a_{jk}, a_{kk} }[/math], которое можно производить параллельно. Вычисление элементов [math]\displaystyle{ a_{jj}, a_{jk}, a_{kk} }[/math] требует порядка [math]\displaystyle{ O(logN) }[/math] операций. Вычисление элементов матриц [math]\displaystyle{ GR, JR }[/math] в параллельной реализации требует [math]\displaystyle{ O(1) }[/math] операций.
Таким образом, при классификации по высоте ЯПФ метод Якоби относится к алгоритмам со сложностью [math]\displaystyle{ O(n^2) }[/math], а при классификации по ширине ЯПФ — к алгоритмам со сложностью [math]\displaystyle{ O(n) }[/math].
Входные и выходные данные алгоритма
Входные данные: плотная матрица [math]\displaystyle{ G }[/math] (элементы [math]\displaystyle{ g_{ij} }[/math]). Дополнительные ограничения:
- [math]\displaystyle{ A= G^TG }[/math] – симметрическая матрица, т. е. [math]\displaystyle{ a_{ij}= a_{ji}, i, j = 1, \ldots, n }[/math].
Объём входных данных: [math]\displaystyle{ n^2 }[/math]
Выходные данные: матрица [math]\displaystyle{ \Sigma = diag(\sigma_{i}) }[/math], где [math]\displaystyle{ \sigma_{i} }[/math] - сингулярные числа, матрица [math]\displaystyle{ U }[/math] левых сингулярных векторов и матрица [math]\displaystyle{ V }[/math] правых сингулярных векторов.
Объём выходных данных: [math]\displaystyle{ 2n^2 + n }[/math] (так как матрица сингулярных чисел [math]\displaystyle{ \Sigma }[/math] - диагональная, то достаточно хранить вектор сингулярных чисел (элементы [math]\displaystyle{ \sigma_i }[/math]), а также две матрицы [math]\displaystyle{ U, V }[/math] левых и правых сингулярных векторов).
Свойства алгоритма
Метод Якоби является самым медленным из имеющихся алгоритмов вычисления сингулярных чисел и сингулярных векторов матрицы. Тем не менее, для некоторых некоторых типов матриц [math]\displaystyle{ G }[/math] он способен вычислять сингулярные числа и векторы намного точнее, чем другие методы. Кроме того, в методе Якоби заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, за счет этого метод достигает высокой относительной точности.
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов для метода Якоби является линейным, вычислительная мощность алгоритма также линейна.
Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности (число итераций зависит от входных данных и порогового значения).
Программная реализация алгоритма
Особенности реализации последовательного алгоритма
Локальность данных и вычислений
Возможные способы и особенности параллельной реализации алгоритма
Масштабируемость алгоритма и его реализации
Масштабируемость алгоритма
Оценки предела масштабируемости алгоритма были рассмотрены в разделе Ресурс параллелизма алгоритма.
Масштабируемость реализации алгоритма
Проведём исследование масштабируемости параллельной реализации метода Якоби нахождения сингулярных чисел и векторов. Исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета.
Набор и границы значений изменяемых реализации алгоритма:
- число процессоров [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
На следующем рисунке приведен график среднего времени работы алгоритма в зависимости от изменяемых параметров запуска.
При увеличении числа процессов эффективность на рассмотренной области изменений параметров запуска падает.
Исследованная параллельная реализация, реализованная в библиотеке Intel MKL
Динамические характеристики и эффективность реализации алгоритма
Выводы для классов архитектур
Существующие реализации алгоритма
Несмотря на то, что метод Якоби является самым медленным из методов вычисления сингулярных чисел и векторов матрицы, он реализован в некоторых пакетах, например, LAPACK, ScaLAPACK, Intel MKL[7].
Литература
- ↑ Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)
- ↑ Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-264)
- ↑ Дж. Деммель «Вычислительная линейная алгебра» (стр. 244-245)
- ↑ Дж. Деммель «Вычислительная линейная алгебра» (стр. 261-263)
- ↑ 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
- ↑ 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
- ↑ Intel® Math Kernel Library, https://software.intel.com/ru-ru/node/521153#CAF8CA08-1C69-48F8-B0C3-B2B27AFA92FA