Преобразования Якоби симметричной матрицы

Метод Якоби всегда является робустным для действительных симметричных матриц. Для матриц размера, скажем больше 10, этот алгоритм более медленный, чем метод QR (см. ниже). Однако алгоритм Якоби значительно проще других более эффективных методов, таким образом он рекомендуется в случаях, когда время выполнения не является критическим.

Базовая матрица ротации Якоби имеет вид:
 

Все диагональные элементы матрицы ротации Якоби равны 1, за исключением двух элементов в строках (и столбцах) p и q. Все внедиагональные элементы равны нулю, за исключением двух элементов, равных s и (-s). Числа c и s являются косинусом и синусом угла поворота f, так что c²+s²=1.

Ротация Якоби преобразует матрицу A согласно формуле:

Для измененных элементов матрицы (с учетом приведенного выше и ее симметрии) получаются явные формулы:

Сходимость метода Якоби можно оценить, рассматривая сумму квадратов внедиагональных элементов S; уравнения преобразований дают для этой суммы после трансформации следующее соотношение:

После цепочки преобразований получается матрица D, диагональная с точностью до машинного нуля. Ее диагональные элементы образуют собственные значения первоначальной матрицы A, поскольку выполняется

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

Для наших целей лучшей стратегией является циклический метод Якоби, где элементы удаляются в строгом порядке. Например, можно просто проходить по строкам: P₁₂ , P₁₃ , ... P_1n ; P₂₃ ,...,P_2n ;... Можно показать, что скорость сходимости будет в общем случае квадратичной, как для оригинального, так и для циклического метода Якоби. Назовем одно такое множество последовательных n(n-1)/2 преобразований Якоби проходом.

Программа, помещенная ниже, основана на алгоритмах из [2, 3]; она содержит еще две тонкости:

• На первых трех проходах мы проводим ротацию относительно индекса (pq) только тогда, когда |a_pq|>e для некоторого порогового значения e = S₀/(5n²), где S₀ -- сумма модулей внедиагональных элементов верхней треугольной матрицы.
• После четырех проходов в том случае, когда |a_pq|<<|a_pp| и |a_pq|<<|a_qq|, полагаем |a_pq|=0 и пропускаем ротацию. Критерием сравнения является превышение диагональных элементов в 10^D+2 раз, где D -- число значащих десятичных цифр.

Типичные матрицы требуют от 6 до 10 проходов для достижения сходимости, или от 3n² до 5n² ротаций. Каждая ротация требует порядка 4n операций умножения и сложения, так что общие затраты имеют порядок от 12n³ до 20n³ операций. При вычислении собственных векторов наряду с собственными значениями, число операций на ротации возрастает с 4n до 6n, что означает превышение только на 50 процентов.

Программа jacobi

#include <math.h>
#include "nrutil.h"    /* Здесь определяются некоторые утилиты типа выделения памяти */

/* Преобразование элементов при ротации */
#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);a[k][l]=h+s*(g-h*tau)
/* максимальное число проходов */
#define MAXSWEEP 50

/* Программа jacobi вычисляет все собственные значения и собственные векторы
   действительной симметричной матрицы a[1...n][1...n]. На выходе разрушаются
   все наддиагональные элементы a. d[1...n] возвращает собственные значения
   матрицы, v[1...n][1...n] -- в столбцах нормализованные собственные векторы,
   nrot -- число ротаций Якоби, потребовавшихся для данной программы.
*/
void jacobi(float **a, int n, float *d, float **v, int *nrot) {
  int j, iq, ip, i;
  float tresh, theta, tau, t, sm, s, h, g, c, *b, *z;
  /* выделить память для временно используемых векторов (декларации в nrutil.h) */
  b=vector(1,n); z=vector(1,n);
  /* инициализировать v как единичную матрицу */
  for(ip=1;ip<=n;ip++) {
    for(iq=1;iq<=n;iq++) v[ip][iq]=0.;
    v[ip][ip]=1.;
  }
  /* инициализировать b диагональю a, z -- нулем */
  for(ip=1;ip<=n;ip++) {b[ip]=a[ip][ip]; z[ip]=0.;}
  /* вначале число ротаций нулевое */
  *nrot=0;
  /* делаем не более MAXSWEEP проходов */
  for(i=1;i<=MAXSWEEP;i++) {
    /* вычисляем сумму модулей внедиагональных элементов */
    for(sm=0.,ip=1;ip<=n;ip++) for(iq=ip+1;iq<=n;iq++) sm += fabs(a[ip][iq]);
   /* диагональная матрица -> нормальный выход */
    if(sm==0.) {
      free_vector(z,1,n); free_vector(b,1,n); return;
    }
    /* пороговое значение элемента для ротации */
    tresh=(i<4 ? 0.2*sm/(n*n) : 0.);
    /* проход осуществляется по строкам, в каждой строке по столбцам */
    for(ip=1;ip<=n-1;ip++) for(iq=ip+1;iq<=n;iq++) {
      /* отследить случай малого элемента после 4 проходов */
      g=100.*fabs(a[ip][iq]);
      if(i>4 && (float)fabs(d[ip]+g)==(float)fabs(d[ip])
             && (float)fabs(d[iq]+g)==(float)fabs(d[iq])) a[ip][iq]=0.;
      /* и случай малого элемента на первых 3 проходах
         (обработать только превысившие порог) */
      else if(fabs(a[ip][iq])>tresh) {
        h=d[ip]-d[iq];
       /* вычислить значение t=s/c по формуле корня квадратного уравнения */
        if((float)(fabs(h)+g)==(float)fabs(h)) t=a[ip][iq]/h;
        else {
          theta=0.5*h/a[ip][iq];
          t=1./(fabs(theta)+sqrt(1.+theta*theta));
          if(theta<0.) t = -t;
        }
       /* вычислить c, s, tau, и др. Изменить диагональ. Обнулить (ip,iq)-элемент */
        c=1./sqrt(1+t*t); s=t*c; tau=s/(1.+c); h=t*a[ip][iq];
        z[ip] -= h; z[iq] += h; d[ip] -= h; d[iq] += h;
        a[ip][iq]=0.;
       /* поворот при 1<=j<ip */
        for(j=1;j<=ip-1;j++) {ROTATE(a,j,ip,j,iq);}
       /* поворот при ip<j<iq */
        for(j=ip+1;j<=iq-1;j++) {ROTATE(a,ip,j,iq,j);}
       /* поворот при iq<j<=n */
        for(j=iq+1;j<=n;j++) {ROTATE(a,ip,j,j,iq);}
       /* добавка для матрицы собственных векторов */
        for(j=1;j<=n;j++) {ROTATE(v,j,ip,j,iq);}
       /* приращение счетчика ротаций */
        ++(*nrot);
      }
    }
    /* добавить до диагонали и реинициализировать z */
    for(ip=1;ip<=n;ip++) {
      b[ip] += z[ip]; d[ip]=b[ip]; z[ip]=0.;
    }
  }
  /* если мы здесь, то число ротаций превысило лимит. Функция nerror (выход с
     диагностикой ошибки) описана декларацией в nrutil.h. */
  nerror("Too many iterations in the routine jacobi");
}

Заметим, что программа, представленная выше, предполагает, что компьютер трактует число ниже машинного нуля как нуль. Если это не так, программа должна быть модифицирована.

На выходе этой программы собственные значения расположены не по порядку величины. Если требуется сортировка, то нижеследующая программа переставит выходные данные программы jacobi или других программ поиска собственных векторов и значений (например, расположенных ниже). Метод сортировки -- прямая вставка -- имеет порядок N² (значительно больше, чем Nlog(N) для эффективных методов, но по сравнению с порядком N³ для процедуры jacobi это ускорение было бы несущественно).

Программа eigsrt

/* Программа eigsrt, получая на входе собственные значения d[1...n] и соответствующие
   собственные векторы v[1...n][1...n], сортирует собственные значения в нисходящем
   порядке, делая соответствующие преобразования в массиве векторов. Метод сортировки
   прямая вставка
*/
void eigsrt(float *d, float **v, int n) {
  int k,j,i;
  float p;
  /* проход по всем значениям собств. чисел, кроме последнего */
  for(i=1;i<n;i++) {
    p=d[k=i];
    /* следует ли делать перестановку и куда */
    for(j=i+1;j<=n;j++) if(d[j]>=p) p=d[k=j];
    /* если следует, произвести ее */
    if(k!=i) {
      d[k]=d[i]; d[i]=p;
      for(j=1;j<=n;j++) {
        p=v[j][i]; v[j][i]=v[j][k]; v[j][k]=p;
      }
   }
}