User:GrishinaAnna/Метод Якоби вычисления сингулярных чисел и векторов

Авторы статьи: Гришина Анна (группа 604), Наталия Чудновская (группа 611)

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

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

Метод Якоби позволяет получить разложение произвольной матрицы [math]\displaystyle{ G }[/math] размера [math]\displaystyle{ n }[/math]×[math]\displaystyle{ n }[/math] в виде [math]\displaystyle{ G=USV^* }[/math], где [math]\displaystyle{ U }[/math] и [math]\displaystyle{ V }[/math] — унитарные матрицы, [math]\displaystyle{ S }[/math] — диагональная матрица с вещественными положительными числами на диагонали. Диагональные элементы матрицы называются сингулярными числами матрицы [math]\displaystyle{ G }[/math], а столбцы матриц и левыми и правыми сингулярными векторами соответственно.

В основе метода лежит утверждение о том, что собственными значениями симметричной матрицы [math]\displaystyle{ G^TG }[/math] являются квадраты сингулярных чисел, а соответствующие ортонормированные собственные векторы – правыми сингулярными векторами. Для вычисления собственных значений и собственных векторов симметричной матрицы [math]\displaystyle{ 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{ G^TG }[/math] не станет достаточно близка к диагональной матрице, после чего вычисляются собственные числа матрицы [math]\displaystyle{ G }[/math], полученной на последнем шаге, которые принимаются за сингулярные числа. Правые сингулярные векторы получаются в результате накопления произведений вращений Якоби.

Метод Якоби для вычисления сингулярных значений и векторов является самым медленным из имеющихся, но тем не менее, интерес к нему сохраняется, потому что для некоторых типов матриц [math]\displaystyle{ G }[/math] он способен вычислять числа и сингулярные векторы намного точнее, чем другие методы. Так, например, метод Якоби вычисляет сингулярные числа матрицы [math]\displaystyle{ G }[/math] с высокой точностью, если [math]\displaystyle{ G }[/math] может быть представлена в виде [math]\displaystyle{ G=DX }[/math], где [math]\displaystyle{ D }[/math] – диагональная матрица, а [math]\displaystyle{ X }[/math] – хорошо обусловлена. Выигрыш метода Якоби в этом случае перед другими алгоритмами объясняется следующим образом: в методе Якоби заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, в то время как другие алгоритмы включают в себя приведение матрицы к двухдиагональной форме, из-за чего и теряют все верные разряды во всех сингулярных числах, кроме старшего.

Следует иметь в виду, что метод Якоби эффективен при не очень больших размерах матрицы – максимум до [math]\displaystyle{ n }[/math] ~ [math]\displaystyle{ 100 }[/math]. Дело не только в том, что метод не будет сходиться за конечное время, дело еще и в излишне большом количестве операций, которое приведет к ошибкам округления и потере точности.

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

Исходные данные: матрица [math]\displaystyle{ G }[/math] размера [math]\displaystyle{ n }[/math]×[math]\displaystyle{ n }[/math].

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

Описание используемой в алгоритме процедуры одностороннего вращения входной матрицы [math]\displaystyle{ G }[/math] One-Sided-Jacobi-Rotation [math]\displaystyle{ (G, j, k) }[/math]:

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

Где [math]\displaystyle{ R(j,k,\theta) }[/math] – матрица вращения Якоби, которая выбирается так, чтобы обнулить пару внедиагональных элементов.

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

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

Построим информационный граф отдельно для процедуры вращения, а затем для всего алгоритма с обращениями к процедуре вращения. Для процедуры вращения в качестве узлов графа выберем точки на равномерной сетке [math]\displaystyle{ 2xN }[/math], где [math]\displaystyle{ N }[/math] – размерность матрицы. Для вычисления коэффициентов в матрице поворота требуется вычислить значения [math]\displaystyle{ a_{ii}, a_{ji}, a_{jj} }[/math] матрицы [math]\displaystyle{ G^TG }[/math]. После того, как посчитаны коэффициенты, необходимо пересчитать значения в [math]\displaystyle{ i }[/math] и [math]\displaystyle{ j }[/math] столбцах матрицы. Эти значения зависят от входных значений соответствующих элементов и рассчитанных ранее значений [math]\displaystyle{ \cos (\theta) }[/math] и [math]\displaystyle{ \sin(\theta) }[/math].

Далее рассмотрим информационный граф для всего алгоритма. Как было показано ранее, процедура вращения в плоскости [math]\displaystyle{ (i,j) }[/math] зависит только от столбцов [math]\displaystyle{ (i,j) }[/math] матрицы [math]\displaystyle{ G }[/math] и рассчитывает новые значения только для [math]\displaystyle{ i }[/math] и [math]\displaystyle{ j }[/math] столбцов, причем новое значение вычисляется для каждого элемента в столбце. С учетом сказанного выше будем строить граф вычислений, рассматривая каждый столбец матрицы G как многомерную переменную и введя обозначение [math]\displaystyle{ G(:,i) }[/math] для [math]\displaystyle{ I-го }[/math] столбца. В зависимости от порядка перебора пар элементов возможны различные графы, параллельные по данным в большей и меньшей степени. При простом переборе for i=1:N; for j=1:N; calc(i,j); end Для каждого последующего элемента требуется столбец из предыдущего внутреннего цикла. С другой стороны, как было показано ранее, процедура односторонних вращений зависит и изменяет только i и j столбцы. Таким образом пары (i1,j1) И (i2,j2) такие, что все i1,i2,j1,j2 – различны можно вычислять независимо.File:Graf 3 chudnovskaya.jpg

Одним из вариантов реализации такого подхода является следующий. Будем строить пары (i,j) специального вида: