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

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

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

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

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

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

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

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

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

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

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

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

Неявный метод Якоби математически эквивалентен применению метода Якоби вычисления собственных значений к матрице [math]\displaystyle{ A = G^TG }[/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{ A = (G^TG) }[/math] можно описать следующим образом:

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

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

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{ 
         \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]\displaystyle{ J = JR(j,k,\theta) }[/math]
        end if
    end if

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

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

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

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

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

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

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

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

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

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

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

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

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

Метод Якоби нахождения сингулярных значений и векторов реализован в немногих пакетах, например, LAPACK, Intel MKL. Это не значит, что алгоритм не представляет никакого интереса для изучения. Для некоторых типов матриц метод способен вычислять сингулярные числа и сингулярные векторы намного точнее, чем другие методы. Но в связи с тем, что алгоритм считается медленным, он не включен во многие известные пакеты.

Литература

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