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

Метод Якоби вычисления сингулярных чисел и векторов
Последовательный алгоритм
Последовательная сложность	[math][/math]
Объём входных данных	[math][/math]
Объём выходных данных	[math][/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math][/math]
Ширина ярусно-параллельной формы	[math][/math]

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

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

Сингулярное разложение было первоначально разработано в дифференциальной геометрии при изучении свойств билинейных форм учеными Эудженио Бельтрами и Камилем Жорданом независимо в 1873 и 1874 годах соответственно. Первое доказательство сингулярного разложения для прямоугольных и комплексных матриц было осуществлено математиками Карлом Эскартом и Гэйлом Янгом в 1936 году.

Сингулярное разложение (Singular Values Decomposition, SVD) является удобным методом при работе с матрицами. Cингулярное разложение показывает геометрическую структуру матрицы и позволяет наглядно представить имеющиеся данные. Сингулярное разложение используется при решении самых разных задач — от приближения методом наименьших квадратов и решения систем уравнений до сжатия и распознавания изображений. Используются разные свойства сингулярного разложения, например, способность показывать ранг матрицы и приближать матрицы данного ранга. Так как вычисления ранга матрицы — задача, которая встречается очень часто, то сингулярное разложение является довольно популярным методом.

Сингулярным разложением матрицы [math]G (n \times n)[/math] называется разложение вида [math]G = U \Sigma V^T[/math] , где [math]U, V[/math] - унитарные матрицы [math]n \times n[/math], [math]\Sigma[/math] - диагональная матрица [math]n \times n[/math] с вещественными положительными числами на главной диагонали. Столбцы матриц [math]U, V[/math] называются соответственно левыми и правыми сингулярными векторами, а значения на диагонали матрицы [math]\Sigma[/math] - сингулярными значениями матрицы [math]G[/math].

Одним из методов нахождения сингулярного разложения матрицы является метод Якоби. Метод Якоби был предложен Карлом Густавом Якоби Якоби в 1846 году и представляет собой итерационный алгоритм вычисления собственных значений и собственных векторов симметричной матрицы. Для вычисления собственных значений необходимо неявно применить метод Якоби к симметричной матрице [math]A = G^T G[/math]. На каждом шаге вычисляется вращение Якоби [math]J[/math], с помощью которого матрица [math]G^TG[/math] неявно пересчитывется в [math]J^TG^TGJ[/math]; вращение выбрано так, чтобы пара внедиагональных элементов из [math]G^TG[/math] обратилась в нули в матрице [math]J^TG^TGJ[/math]. При этом ни [math]G^TG[/math], ни [math]J^TG^TGJ[/math] не вычисляются в явном виде; вместо них вычисляется матрица [math]GJ[/math]. Поэтому алгоритм называется методом односторонних вращений.

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

Исходные данные: Матрица [math]G[/math] размера [math](n \times n)[/math] над полем вещественных или комплексных чисел.

Вычисляемые данные:

• [math]\sigma_i[/math] - сингулярные числа матрицы [math]G[/math],
• [math]U[/math] - матрица левых сингулярных векторов,
• [math]V[/math] - матрица правых сингулярных векторов.

Описание алгоритма:

Шаг 1 Вычисляем матрицу [math]G'[/math] применением одностороннего вращения Якоби к исходной матрице [math]G[/math].

Шаг 2 Повторяем Шаг 1 пока матрица [math]G'^TG[/math] не станет достаточно близка к диагональной.

Шаг 3 Вычисляем по полученной матрице [math]G'[/math] сингулярные числа [math]\sigma_i[/math] и сингулярные векторы [math]U, V[/math].

Основные формулы метода:

[math]\sigma_i = \parallel G'(:,i)\parallel_2[/math] (2-я норма [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]a_{ij}[/math] - элемент, расположенный в [math]i[/math]-ой строке и [math]j[/math]-ом столбце матрицы [math]G^TG[/math].

[math]\tau = (a_{jj} - a_{kk})/(2 \cdot a_{jk})[/math]

[math]t = sign(\tau)/(|\tau|+\sqrt{1+\tau^2})[/math]

[math]c = 1/\sqrt{1+t^2}[/math], где [math]c = cos\theta[/math]

[math]s = c \cdot t[/math], где [math]s = sin\theta[/math]

[math]J = J \cdot R(j,k,\theta)[/math]

Здесь [math]R(j,k,\theta)[/math] - матрица вращений Якоби, которая имеет следующий вид:

                                                                                 [math] \begin{matrix} & & & & j & & k & & & \\ \end{matrix} [/math]

[math] R(j,k,\theta) = \begin{matrix} \\ \\ \\ j \\ \\ k \\ \\ \\ \\ \end{matrix} [/math] [math] \begin{pmatrix} & 1 & & & & & & & & \\ & & 1 & & & & & & & \\ & & & \ddots & & & & & & \\ & & & & \cos(\theta) & & -\sin(\theta) & & & \\ & & & & & \ddots & & & & \\ & & & & \sin(\theta) & & \cos(\theta) & & & \\ & & & & & & & \ddots & & \\ & & & & & & & & 1 & \\ & & & & & & & & & 1 \\ \end{pmatrix} [/math]

Одностороннее вращение Якоби в координатной плоскости [math]i, j[/math] вычисляется в том случае, если элемент [math]a_{ij}[/math] удовлетворяет условию: [math]|a_{ij}| \geq \epsilon \sqrt{a_{ii}a_{jj}}[/math], где [math]\epsilon[/math] -

Критерий останова:

Для проверки близости матрицы [math]G'^TG[/math] к диагональной используется следующее неравенство:

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

Вычислительное ядро метода Якоби вычисления сингулярных чисел и векторов можно составить из множественных вызовов процедуры одностороннего вращения Якоби One-Sided-Jacobi-Rotation(G,j,k).

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

• Вычисление элементов [math]a_{jj}, a_{jk}, a_{kk}[/math], где [math]a_{ij} = (G^TG)_{ij}[/math].
• Вычисление матрицы вращения Якоби [math]R(j,k,\theta)[/math].
• Применение вращения Якоби к матрицам [math]G[/math] и [math]J[/math].

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

Алгоритм вычисления сингулярных чисел и векторов Якоби состоит из множественных обращений к процедуре одностороннего вращения Якоби. На каждой итерации цикла происходит вызов процедуры One-Sided-Jacobi-Rotation(G, j, k) для всех [math]j \in [1, n-1][/math] и [math]k \in [j+1, n][/math]. Таким образом, на каждой итерации цикла процедура One-Sided-Jacobi-Rotation(G, j, k) вызывается для всех наддиагональных элементом матрицы [math]G^TG[/math], то есть всего [math]\frac{n(n-1)}{2}[/math] раз. Количество итераций цикла непосредственно зависит от матрицы [math]G[/math], выход из цикла происходит по достижении критерия останова (описан в разделе "Математическое описание алгоритма".)

1.5 Схема реализации последовательного алгоритма

Схему реализации последовательного алгоритма можно описать следующиv псевдокодом:

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(G, j, k)
    end for
  end for
пока [math]G^TG[/math] не станет достаточно близка к диагональной матрице.
Положить [math]\sigma_i = \parallel G'(:,i)\parallel_2[/math] (2-я норма [math]i[/math]-го столбца в [math]G'[/math]).
Положить [math]U = [u_1, \dots, u_n][/math], где [math] u_i = G'(:,i)/\sigma_i[/math].
Положить [math]V = J[/math], где [math]J[/math] - накопленное произведение вращений Якоби.

Процедура One-Sided-Jacobi-Rotation(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]\tau = {a_{jj} - a_{kk}}/{2a_{jk}}[/math]
        [math]t = {sign(\tau)}/(|\tau| + \sqrt{1 + \tau^2})[/math]
        [math]c = 1/(\sqrt{1 + t^2})[/math]
        [math]s = ct[/math]
        [math]G = GR(j,k,\theta) \dots [/math]где [math]c=cos\theta, s=sin\theta[/math]
        if нужны правые сингулярные векторы 
            [math]J = JR(j,k,\theta)[/math]
        end if
    end if

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

Одним из преимуществ алгоритма Якоби вычисления сингулярных чисел и векторов перед другими алгоритмами является его параллелизм.

Алгоритм Якоби позволяет производить параллельно следующие вычисления:

* Вычисление элементов матрицы [math]G^TG[/math]: необходимо вычислить три элемента матрицы [math]G^TG[/math]: [math]a_{jj}, a_{jk}, a_{kk}[/math]. Вычисление одного элемента не зависит от вычисления других, поэтому их можно производить параллельно.

Сложность операции: сложность вычисления одного элемента сравнима с [math]O(n)[/math], при условии параллельного вычисления всех трех элементов суммарная сложность операции остается сравнимой с [math]O(n)[/math].

* Вычисление элементов матриц [math]GR[/math], [math]JR[/math]: умножение матрицы [math]G[/math] на матрицу вращений Якоби ([math]R(j,k\theta)[/math]) изменяет всего два столбца исходной матрицы ([math]j[/math]-й и [math]k[/math]-й). В сумме необходимо вычислить [math]4n[/math] элемента ([math]2n[/math] элементов матрицы [math]G[/math] и [math]2n[/math] элементов матрицы [math]J[/math]), вычисление всех элементов можно производить параллельно.

Сложность операции: сложность вычисления одного элемента сравнима с [math]O(n)[/math], при условии параллельного вычисления всех элементов суммарная сложность операции становится сравнимой с [math]O(n^2)[/math].

* Вычисление вращений Якоби для всех наддиагональных элементов: как описано выше, одно вращения Якоби изменяет всего два столбца исходной матрицы, поэтому параллельно могут быть вычислены вращения для [math]\lceil n/2 \rceil[/math] столбцов матрицы.

Сложность операции: Вычисление вращения Якоби для всех наддиагональных элементов состоит из [math]\frac{n(n-1)}{2}[/math] итераций. При выборе столбцов [math]j[/math] и [math]k[/math] таким образом, чтобы параллельно можно было выполнить [math]\lceil n/2 \rceil[/math] вращений Якоби (вариант выбора итераций описан в разделе Информационный граф), количество последовательных итераций сокращается до [math]n-1[/math]. Суммарная сложность операции при условии параллельного вычисления становится сравнимой с [math]O(n^3)[/math].

3 Литература

<references \>