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

Метод Якоби вычисления сингулярных чисел и векторов
Последовательный алгоритм
Последовательная сложность	[math]O(n^3)[/math]
Объём входных данных	[math]n^2[/math]
Объём выходных данных	[math]2n^2+n[/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]G = G \cdot R(j,k,\theta)[/math]

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

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

                                                                [math]j[/math]                          [math]k[/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 \varepsilon \sqrt{a_{ii}a_{jj}}[/math].

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

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

[math]off(A) = \sqrt{\sum_{1\leq j\lt k \leq n} a^2_{jk}} \leq \varepsilon[/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)  [/math], где [math]c=cos\theta, s=sin\theta[/math] 
        if нужны правые сингулярные векторы 
            [math]J = JR(j,k,\theta)[/math]
        end if
    end if

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

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

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

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

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

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

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

Опишем граф алгоритма на макро- и микро-уровне.

При одном вращении Якоби изменяется всего 2 столбца матрицы [math]G[/math]. Таким образом, есть возможность одновременно выполнять [math]\lceil \frac{n}{2} \rceil[/math] вращений.

Рисунок 1. Информационный граф на макроуровне.

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

• Вычисление элементов матрицы [math]G^TG[/math]: необходимо вычислить три элемента матрицы [math]G^TG[/math]: [math]a_{jj}, a_{jk}, a_{kk}[/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]), вычисление всех элементов можно производить параллельно.

Рисунок 2 иллюстрирует ЯПФ для процедуры одностороннего вращения Якоби. Каждый серый овальный прямоугольник объединяет вычисления на одном уровне ЯПФ, т.е. вычисления, которые могут быть произведенные параллельно. Каждый розовый прямоугольник представляет собой набор базовых математических операций типа сложение, умножение, деление и т.д.

Рисунок 2. Информационный граф процедуры одностороннего вращения Якоби.

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].

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

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

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

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

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

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

3 Литература

<references \>