Участник:Galkina/Метод Якоби вычисления собственных значений симметричной матрицы

Основные авторы описания: А.С.Галкина, И.А.Плахов

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

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

Метод Якоби — итерационный алгоритм для вычисления собственных значений и собственных векторов вещественной симметричной матрицы. Карл Густав Якоб Якоби, в честь которого назван этот метод, предложил его в 1846 году, хотя использоваться он начал только в 1950-х годах с появлением компьютеров. Суть алгоритма заключается в том, чтобы по заданной симметрической матрице [math]A = A_0[/math] построить последовательность ортогонально подобных матриц [math]A_1,A_2,...[/math]. Эта последовательность сходится к диагональной матрице, на диагонали которой стоят собственные значения.

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

Исходные данные: симметрическая матрица [math]A[/math] (элементы [math]a_{ij}[/math]).

Вычисляемые данные: диагональная матрица [math]\Lambda[/math] (элементы [math]\lambda_{ij}[/math]).

Матрица [math]A_{i+1}[/math] получается из [math]A_i[/math] по формуле [math]A_{i+1}={J_i}^TA_iJ_i[/math], где [math]J_i[/math] — ортогональная матрица, называемая вращением Якоби. При подходящем выборе [math]J_i[/math] матрица [math]A_m[/math] для больших m будет близка к диагональной матрице [math]\Lambda[/math].

Матрица [math]J_i[/math] выбирается так, чтобы сделать нулями пару внедиагональных элементов матрицы [math]A_{i+1}[/math].

                                                 [math]j[/math]                         [math]k[/math]

[math] J_i = \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]A_{i+1}[/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]

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

[math] \begin{bmatrix} a_{jj}^{(i+1)} & a_{jk}^{(i+1)} \\ a_{kj}^{(i+1)} & a_{kk}^{(i+1)} \end{bmatrix} = \begin{bmatrix} \cos(\theta) & -\sin(\theta) \\ \sin(\theta) & \cos(\theta) \end{bmatrix}^T \begin{bmatrix} a_{jj}^{(i)} & a_{jk}^{(i)} \\ a_{kj}^{(i)} & a_{kk}^{(i)} \end{bmatrix} \begin{bmatrix} \cos(\theta) & -\sin(\theta) \\ \sin(\theta) & \cos(\theta) \end{bmatrix} = \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix}[/math]

где [math]\lambda_1[/math] и [math]\lambda_2[/math] - собственные значения подматрицы

[math] \begin{bmatrix} a_{jj}^{(i+1)} & a_{jk}^{(i+1)} \\ a_{kj}^{(i+1)} & a_{kk}^{(i+1)} \end{bmatrix} [/math]

Опуская индекс [math](i)[/math], с учетом симметрии получим:

[math] \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix} = \begin{bmatrix} a_{jj}c^2 + a_{kk}s^2 + 2sca_{jk} & sc(a_{kk}-a_{jj})+a_{jk}(c^2-s^2) \\ sc(a_{kk}-a_{jj})+a_{jk}(c^2-s^2) & a_{jj}s^2 + a_{kk}c^2 - 2sca_{jk} \end{bmatrix} [/math]

Приравнивая внедиагональный элемент нулю, получаем

[math] \frac{a_{jj}^{(i)} - a_{kk}^{(i)}}{2 a_{jk}^{(i)}} = \frac{c^2 - s^2}{2sc} = \frac{\cos 2\theta}{\sin 2\theta} = \operatorname{tg}(2\theta) \equiv \tau [/math].

Положим [math]t = \frac{s}{c} = \operatorname{ctg}(\theta)[/math] и заметим, что [math]t^2 - 2t\tau + 1 = 0[/math]. Решая квадратное уравнение, находим [math]t = \frac{\operatorname{sign}(\tau)}{|\tau| + \sqrt{1+\tau^2}}, c = \frac{1}{\sqrt{1+t^2}}, s = tc[/math].

Выбор параметров [math]j[/math] и [math]k[/math] производится путем построчного циклического обхода внедиагональных элементов матрицы [math]A[/math].

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

Вычислительное ядро составляют множественные вычисления элементов матрицы [math]a_{jm}^{(i+1)} = a_{mj}^{(i+1)}[/math]   и   [math]a_{km}^{(i+1)} = a_{mk}^{(i+1)}[/math]   в процессе применения матрицы поворота [math]J[/math] к матрице [math]A[/math]:

[math]\begin{align} 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 \\ \end{align}[/math]

каждое из которых повторяется по [math] (n-2) [/math] раза.

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

Основную часть метода составляет процедура применения вращения к матрице [math]A[/math], которая в дальнейшем будет обозначена как Jakobi-Rotation[math](A,j,k)[/math]

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

Алгоритм можно описать следующим образом:

repeat
    for [math]j=1[/math] to [math]n-1[/math]
        for [math]k=j+1[/math] to [math]n[/math]
            выполнить процедуру Jakobi-Rotation[math](A,j,k)[/math]
        end for
    end for
пока [math]A[/math] не достаточно близка к диагональной матрице

Процедура Jakobi-Rotation[math](A,j,k)[/math] — это следующий алгоритм:

if [math]|a_{jk}|[/math] не слишком мал
    [math]\begin{align}
      \tau &= \frac{a_{jj}-a_{kk}}{2\,a_{jk}} \\
      t &= \frac{sign(\tau)}{|\tau|+\sqrt{1+\tau^2}} \\
      c &= \frac{1}{\sqrt{1+t^2}} \\
      s &= c\,t \\
      A &= R^T(j,k,\theta)\cdot A\cdot R(j,k,\theta), \qquad J_i = R(j,k,\theta),\ c = \cos \theta,\ s = \sin \theta 
     \end{align}[/math]
end if

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

Для вычисления собственных значений вещественной симметричной матрицы порядка n методом Якоби требуется:

Для вычислительного ядра —

• [math] 2\,n(n-1)(n-2) [/math] умножений,
• [math] n(n-1)(n-2) [/math] сложений.

Для остальной части алгоритма —

• [math] \frac{3\,n(n-1)}{2} [/math] умножений,
• [math] \frac{3\,n(n-1)}{2} [/math] делений,
• [math] \frac{5\,n(n-1)}{2} [/math] сложений (вычитаний),
• [math] n(n-1) [/math] операций извлечения квадратного корня.

Умножения и сложения (вычитания) составляют основную часть алгоритма.

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

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

Входные данные: матрица [math]A[/math] (элементы [math]a_{ij}[/math]). Дополнительные ограничения:

• [math]A[/math] – симметрическая матрица, т. е. [math]a_{ij}= a_{ji}, i, j = 1, \ldots, n[/math].

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

Выходные данные: диагональная матрица [math]\Lambda[/math] (элементы [math]\lambda_{ii}[/math]).

Объём выходных данных: [math] n [/math] (достаточно хранить только диагональные элементы).