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

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

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

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

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

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

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

Одним из методов нахождения сингулярного разложения матрицы является метод Якоби. Метод Якоби был предложен Карлом Густавом Якоби Якоби в 1846 году и представляет собой итерационный алгоритм вычисления собственных значений и собственных векторов симметричной матрицы. Для вычисления собственных значений необходимо неявно применить метод Якоби к симметричной матрице [math]\displaystyle{ A = G^T G }[/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{ (n \times n) }[/math] над полем вещественных или комплексных чисел.

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

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

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

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

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

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

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

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

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

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

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

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

[math]\displaystyle{ G = G \cdot R(j,k,\theta) }[/math]

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

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

                                                                [math]\displaystyle{ j }[/math]                          [math]\displaystyle{ k }[/math]            

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

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

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

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

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

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

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

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

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

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

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

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

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

Процедура One-Sided-Jacobi-Rotation(G, j, k):

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

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

Как описано выше, для вычисления сингулярных чисел и векторов методом Якоби необходимо [math]\displaystyle{ \frac{n(n-1)}{2} }[/math] выполнение процедуры одностороннего вращения Якоби.

Cложность процедуры One-Sided-Jacobi-Rotation(G, j, k):

1. Вычисление элементов [math]\displaystyle{ a_{jj}, a_{jk}, a_{kk} }[/math], где [math]\displaystyle{ a_{ij} = (G^TG)_{ij} }[/math]. Для вычисления одного элемента [math]\displaystyle{ a_{ij} }[/math] достаточно перемножить только векторы [math]\displaystyle{ G^T(i,:) }[/math] и [math]\displaystyle{ G(:,j) }[/math]. Соответственно потребуется [math]\displaystyle{ 3(2n-1) }[/math] операций умножения и сложения.
2. Вычисление значений [math]\displaystyle{ \tau, t, c, s }[/math] для матрицы вращения Якоби [math]\displaystyle{ R(j,k,\theta) }[/math] требует порядка [math]\displaystyle{ O(1) }[/math] операций.
3. Вычисление [math]\displaystyle{ G = GR(j,k,\theta) }[/math] и [math]\displaystyle{ J = JR(j,k,\theta) }[/math]. Можно заметить, что [math]\displaystyle{ G_{il} }[/math] при [math]\displaystyle{ l \ne j,k }[/math] при умножении на матрицу [math]\displaystyle{ R(j,k,\theta) }[/math] не изменяются, и необходимо произвести вычисление только столбцов [math]\displaystyle{ j }[/math] и [math]\displaystyle{ k }[/math]. При этом для вычисления одного элемента требуется 3 операции, так как соответствующий столбец матрицы [math]\displaystyle{ R(j,k,\theta) }[/math] имеет только 2 ненулевых элемента. Следовательно, для вычисления [math]\displaystyle{ G = GR(j,k,\theta) }[/math] необходимо выполнить [math]\displaystyle{ 6n }[/math] операций и столько же для [math]\displaystyle{ J = JR(j,k,\theta) }[/math].

Таким образом, одно вращение Якоби требует порядка [math]\displaystyle{ O(n) }[/math] операций.

Последовательная сложность алгоритма Якоба равна [math]\displaystyle{ \frac{n(n-1)}{2} \cdot O(n) = O(n^3) }[/math].

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

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

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

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

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

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

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

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

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

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

Входные данные: Матрица [math]\displaystyle{ G }[/math] размера [math]\displaystyle{ (n \times n) }[/math].

Объём входных данных: [math]\displaystyle{ n^2 }[/math]

Выходные данные:

• [math]\displaystyle{ \sigma_i }[/math] - сингулярные числа матрицы [math]\displaystyle{ G }[/math] ([math]\displaystyle{ i = 1,\dots,n }[/math]),
• [math]\displaystyle{ U }[/math] - матрица размера [math]\displaystyle{ (n \times n) }[/math] левых сингулярных векторов,
• [math]\displaystyle{ V }[/math] - матрица размера [math]\displaystyle{ (n \times n) }[/math] правых сингулярных векторов.

Объём выходных данных: [math]\displaystyle{ 2n^2+n }[/math]

Литература

<references \>