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

Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Frolov и ASA.

Метод Якоби вычисления сингулярных чисел и векторов
Последовательный алгоритм
Последовательная сложность	$O(n^3)$
Объём входных данных	$n^2$
Объём выходных данных	$2n^2 + n$
Параллельный алгоритм
Высота ярусно-параллельной формы	$O(n^2)$
Ширина ярусно-параллельной формы	$O(n)$

Авторы статьи: Почернина Елена, группа 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)

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

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

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

Неявный метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице $A = G^TG$. Это значит, что на каждом шаге вычисляется вращение Якоби $J$, с помощью которого матрица $G^TG$ неявно пересчитывется в $J^TG^TGJ$; вращение выбрано так, чтобы пара внедиагональных элементов из $G^TG$ обратилась в нули в матрице $J^TG^TGJ$. При этом ни $G^TG$, ни $J^TG^TGJ$ не вычисляются в явном виде; вместо них вычисляется матрица $GJ$. Поэтому алгоритм называется методом односторонних вращений.

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

1.2 Математическое описание алгоритма

Исходные данные: матрица $G$ (элементы $g_{ij}, i, j = 1, \ldots, n$).

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

Основной алгоритм заключается в следующем:

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

$\sigma_i = ||G(:,i)||_2$;

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

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

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

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

Здесь $R(j,k,\theta)$ - матрица вращений Якоби. Матрица вращений выбирается так, чтобы обнулить пару внедиагональных элементов.

                                                                                 $\begin{matrix} & & & & j & & k & & & \\ \end{matrix}$

$J_i = R(j,k,\theta) = \begin{matrix} \\ \\ \\ j \\ \\ k \\ \\ \\ \\ \end{matrix}$ $\begin{bmatrix} & 1 & & & & & & & & \\ & & 1 & & & & & & & \\ & & & \ddots & & & & & & \\ & & & & \cos(\theta) & & -\sin(\theta) & & & \\ & & & & & \ddots & & & & \\ & & & & \sin(\theta) & & \cos(\theta) & & & \\ & & & & & & & \ddots & & \\ & & & & & & & & 1 & \\ & & & & & & & & & 1 \\ \end{bmatrix}$

При совпадении диагональных элементов $(a_{jj} = a_{kk})$ считаем $\theta = \frac{\pi}{4}$.

Для завершения описания алгоритма, нужно указать критерий останова, то есть объяснить смысл условия "если $a_{jk}$ не слишком мал". Подходящим критерием останова является условие $\begin{align} |a_{jk}| \ge \varepsilon\sqrt{a_{jj}a_{kk}} \end{align}$.

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

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

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

$\sigma_i = ||G(:,i)||_2$

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

$u_i = G(:,i)/\sigma_i$

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

$V = J, J = JR(j, k, \theta)$

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

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

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

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

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

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

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

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

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

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

Метод Якоби для нахождения сингулярных чисел и векторов матрицы $G$ также использует метод односторонних вращений (процедура 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

1.6 Последовательная сложность алгоритма

Для вычисления сингулярных чисел и векторов матрицы порядка $n$ в последовательном варианте на каждой итерации алгоритма вызывается процедура One-Sided-Jacobi-Rotation$(G,j,k)$, сложность которой равна $O(n)$. Так как выбор индексов $j$ и $k$ осуществляется путем перебора внедиагональных элементов матрицы, для полного прохода потребуется $\frac{n(n-1)}{2}$ итераций.

Таким образом, метод Якоби вычисления сингулярных чисел и векторов при классификации по последовательной сложности относится к алгоритмам с кубической сложностью.

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

Граф алгоритма состоит из трех групп вершин:

Первая группа вершин расположена в двумерной области, соответствующая ей операция - вычисление $a_{jj}, a_{jk}, a_{kk}$. Естественно введённые координаты области таковы:

• $j$ — меняется в диапазоне от $1$ до $n-1$, принимая все целочисленные значения;
• $k$ — меняется в диапазоне от $j+1$ до $n$, принимая все целочисленные значения.

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

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

• 

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

• 

  Рисунок 2. Граф вычисления значений элементов $a_{jj}, a_{jk}, a_{kk}$.

• 

  Рисунок 3. Граф одной итерации алгоритма с вычислением матрицы вращения $R$ отображением входных и выходных данных.

• 

  Рисунок 4. Граф вычисления значений элементов матрицы $GR$

• 

  Рисунок 5. Граф вычисления значений элементов матрицы $JR$

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

В данном алгоритме существует возможность параллельного вычисления вращений Якоби. Максимальное количество вращений в данном методе равно $\frac{n(n-1)}{2}$. В параллельной реализации, если $n$ - четное, можно выполнить $\frac{n}{2}$ вращений параллельно, если $n$ - нечетное, то $\frac{n-1}{2}$ вращений.

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

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

1.9 Входные и выходные данные алгоритма

Входные данные: плотная матрица $G$ (элементы $g_{ij}$). Дополнительные ограничения:

• $A= G^TG$ – симметрическая матрица, т. е. $a_{ij}= a_{ji}, i, j = 1, \ldots, n$.

Объём входных данных: $n^2$

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

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

1.10 Свойства алгоритма

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

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

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

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

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

2.4.1 Масштабируемость алгоритма

Оценки предела масштабируемости алгоритма были рассмотрены в разделе Ресурс параллелизма алгоритма.

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

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

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

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

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

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

При запуске использовалась версия 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

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

Рисунок 7. Параллельная реализация метода Якоби вычисления сингулярных чисел и векторов. Изменение производительности в зависимости от числа процессов и размера матрицы

Рисунок 8. Параллельная реализация метода Якоби вычисления сингулярных чисел и векторов. Изменение эффективности в зависимости от числа процессов и размера матрицы

Построим оценки масштабируемости реализации метода Якоби вычисления сингулярных чисел и векторов:

• По числу потоков: при увеличении числа потоков эффективность на рассмотренной области увеличивается при переходе от двух потоков к трем. Дальнейшее увеличение количества потоков приводит к снижению эффективности.
• По размеру задачи: наибольшая производительность достигается на матрицах размера 300x300 и 400x400. Дальнейшее увеличение размера матрицы не приводит к увеличению производительности.

Исследованная параллельная реализация, реализованная в библиотеке Intel MKL

2.7 Существующие реализации алгоритма

Несмотря на то, что метод Якоби является самым медленным из методов вычисления сингулярных чисел и векторов матрицы, он реализован в некоторых пакетах, например, LAPACK, ScaLAPACK, Intel MKL.