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

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

Авторы статьи: Давлетшина Александра (группа 619), Зайцева Александра (группа 601)

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

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

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

• левые сингулярные векторы матрицы А — это собственные векторы матрицы [math]AA^*[/math];

• правые сингулярные векторы матрицы A — это собственные векторы матрицы [math]A^*A[/math].

Геометрическая интерпретация сингулярного разложения заключается в следующем факте из геометрии: Образом любого линейного преобразования, заданного с помощью матрицы [math]m[/math]x[math]n[/math], примененного к единичной сфере является гиперэллипсоид.

Из многочисленных разложений матриц, наиболее удобным является сингулярное разложение матрицы в виде: [math]A=U\Sigma V^*[/math], где [math]U, V[/math] – унитарные матрицы, а [math]\Sigma[/math] – диагональная матрица с вещественными положительными числами на диагонали .

Диагональные элементы матрицы [math]\Sigma[/math] называются сингулярными числами матрицы [math]A[/math], а столбцы матриц [math]U,V[/math] левыми и правыми сингулярными векторами соответственно.

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

Неявный метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице [math]A = G^TG[/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]g_{ij}, i, j = 1, \ldots, n[/math]).

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

симметричную матрицу [math]A = (G^TG)[/math] можно привести к диагональному виду: [math]RG^TGR^T[/math], где

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

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

Обозначим [math]s = \sin \theta[/math] и [math]c = \cos \theta[/math].

Перемножив матрицы [math]RG^TGR^T[/math], получим следующие выражения:

[math]\begin{align} a_{jj}^{(i+1)} &= c^2\, a_{jj}^{(i)} - 2\, s c \,a_{jk}^{(i)} + s^2\, a_{kk}^{(i)} \\ a_{kk}^{(i+1)} &= s^2 \,a_{jj}^{(i)} + 2 s c\, a_{jk}^{(i)} + c^2 \, a_{kk}^{(i)} \\ a_{jk}^{(i+1)} &= a_{kj}^{(i+1)} = (c^2 - s^2 ) \, a_{jk}^{(i)} + s c \, (a_{kk}^{(i)} - a_{jj}^{(i)} ) \\ a_{jm}^{(i+1)} &= a_{mj}^{(i+1)} = c \, a_{jm}^{(i)} - s \, a_{km}^{(i)} & m \ne j,k \\ a_{km}^{(i+1)} &= a_{mk}^{(i+1)} = s \, a_{jm}^{(i)} + c \, a_{km}^{(i)} & m \ne j,k \\ a_{ml}^{(i+1)} &= a_{ml}^{(i)} &m,l \ne j,k \end{align}[/math]

Т.к внедиагональные элементы должны быть равны 0, то выбираем [math]\theta[/math] так, чтобы [math]a_{jk}^{(i+1)} = 0[/math] и [math]a_{kj}^{(i+1)} = 0[/math].

Тогда решая уравнение получаем, что

[math] \begin{align} \tau &= {a_{jj} - a_{kk}}/{2a_{jk}} \\ t &= {sign(\tau)}/(|\tau| + \sqrt{1 + \tau^2}) \\ c &= 1/(\sqrt{1 + t^2}) \\ s &= ct \\ \end{align}[/math]

Остановка вычислений алгоритма наступает, когда внедиагональные элементы становятся меньше некоторого наперед заданного числа [math]\varepsilon[/math].

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

Вычислительное ядро алгоритма составляет непосредственно вычисление одностороннего вращения Якоби, т.е. вычисление элементов матрицы [math]g_{jm}^{(i+1)} = g_{mj}^{(i+1)}[/math] и [math]g_{km}^{(i+1)} = g_{mk}^{(i+1)}[/math], [math]m \ne j,k[/math], а также вычисления элементов [math]g_{jj}^{(i+1)} [/math] и [math]g_{kk}^{(i+1)} [/math], получившихся после умножения промежуточной матрицы на матрицу вращений Якоби:

[math]\begin{align} g_{jj}^{(i+1)} &= c^2\, a_{jj}^{(i)} - 2\, s c \,g_{jk}^{(i)} + s^2\, g_{kk}^{(i)} \\ g_{kk}^{(i+1)} &= s^2 \,a_{jj}^{(i)} + 2 s c\, g_{jk}^{(i)} + c^2 \, g_{kk}^{(i)} \\ g_{jk}^{(i+1)} &= g_{kj}^{(i+1)} = (c^2 - s^2 ) \, g_{jk}^{(i)} + s c \, (g_{kk}^{(i)} - g_{jj}^{(i)} ) \\ g_{jm}^{(i+1)} &= g_{mj}^{(i+1)} = c \, g_{jm}^{(i)} - s \, g_{km}^{(i)} & m \ne j,k \\ g_{km}^{(i+1)} &= g_{mk}^{(i+1)} = s \, g_{jm}^{(i)} + c \, g_{km}^{(i)} & m \ne j,k \\ g_{ml}^{(i+1)} &= g_{ml}^{(i)} &m,l \ne j,k \end{align}[/math]

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

Основную часть метода составляет процедура односторонних вращений матрицы, которая в свою очередь состоит из:

1. Вычисление [math]\tau, c, s, t[/math].
2. Одностороннее вращение матрицы.

Что требует [math]O(n)[/math] операций.

После этого необходимо:

1. Вычислить сингулярные числа [math]\sigma_i[/math]. Что требует [math]O(n)[/math] операций.
2. Вычислить [math]U[/math] и [math]V[/math]. Что требует [math]O(n^2)[/math] операций.

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

Применение метода Якоби к матрице [math]A = (G^TG)[/math] можно описать следующим образом:

Если [math]G = U\Sigma V^T[/math] и [math]\Sigma = diag(\sigma_{i})[/math], то алгоритм вычисляет сингулярные числа [math]\sigma_{i}[/math], матрицу [math]U[/math] левых сингулярных векторов и матрицу [math]V[/math] правых сингулярных векторов по следующей схеме:

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

Где процедура 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 &= {a_{jj} - a_{kk}}/{2a_{jk}} \\
         t &= {sign(\tau)}/(|\tau| + \sqrt{1 + \tau^2}) \\
         c &= 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.8 Ресурс параллелизма алгоритма

Для выполнения одной итерации метода Якоби требуется последовательно выполнить несколько ярусов:

1. вычислить [math]\tau[/math]
2. вычислить [math]t[/math]
3. вычислить [math]c[/math]
4. выполнить поворот матрицы [math]G[/math] с параллельным выполнением [math]2(n-2)[/math] операций вычисления элементов матрицы [math]g_{jm}^{(i+1)} = g_{mj}^{(i+1)}[/math] и [math]g_{km}^{(i+1)} = g_{mk}^{(i+1)}[/math], [math]m \ne j,k[/math], а также вычисления элементов [math]g_{jj}^{(i+1)} [/math] и [math]g_{kk}^{(i+1)} [/math].

Суммарное количество итераций равно [math]\frac{n(n-1)}{2}[/math] (Различные варианты выбора [math]j[/math] и [math]k[/math]).

Асимптотически метод сходится квадратично.

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

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

Входные данные: плотная матрица [math]G[/math] размера [math]n[/math]х[math]n[/math], где [math]A= G^TG[/math] – симметрическая матрица.

Объём входных данных: [math]n^2[/math] (т.к. для хранения матрицы [math]G[/math] используем двумерный массив размера [math]n[/math]х[math]n[/math])

Выходные данные: [math]\sigma_{i}[/math] - сингулярные числа, матрица [math]U[/math] левых сингулярных векторов и матрица [math]V[/math] правых сингулярных векторов.

Объём выходных данных: [math]2n^2 + n[/math] (т.к. необходимо хранить вектор сингулярных чисел длины [math]n[/math], а также две матрицы [math]U, V[/math] левых и правых сингулярных векторов размера [math]n[/math]х[math]n[/math]).

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

Метод Якоби является самым медленным из имеющихся алгоритмов вычисления сингулярных чисел и сингулярных векторов матрицы. Тем не менее, для матриц [math]G[/math], которые могут быть представлены в виде [math]G=DX[/math], где [math]D[/math] - диагональная, а [math]X[/math] - хорошо обусловлена, он способен вычислять сингулярные числа и вектора намного точнее, чем другие методы.

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

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

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

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

Метод Якоби нахождения сингулярных значений и векторов реализован в таких пакетах, как LAPACK, LINPACK/EISPACK, Intel MKL .

Так же функция Якоби реализована в пакете Matlab.

Т. к. алгоритм считается медленным, он не включен во многие известные пакеты.

3 Литература

1. Дж. Деммель Вычислительная линейная алгебра. Изд. Мир, 2001.