Собственные вектора и значения трехдиагональной матрицы

Программы, основанные на этих идеях, имеются в [2, 3]. Однако, если требуется только небольшая часть от всей совокупности собственных векторов и значений, значительно более эффективным будет метод факторизации, помещенный ниже.

Алгоритмы QR и QL

Теперь рассмотрим матрицу, сформированную произведением этих факторов в обратном порядке:

Можно проверить, что QR-преобразование сохраняет следующие свойства матрицы: симметрию, трехдиагональность и форму Гессенберга (определенную ниже).

Не обязательно выбирать один из факторов A в виде верхней треугольной матрицы: можно точно также выбрать нижнюю треугольную. Это приведет к QL-алгоритму, поскольку разложение A будет записано в форме

Вспомним, что редукцию Хаусхолдера мы проводили, начиная с n-ого (последнего) столбца исходной матрицы. Для уменьшения ошибки округления мы советовали расположить ее наибольшие элементы ближе к правому нижнему углу. Если мы после этого хотим диагонализировать получившуюся трехдиагональную матрицу, то алгоритм QL приведет к меньшей ошибке округления, чем QR, таким образом далее мы будем говорить об алгоритме QL.

Алгоритм QL состоит из цепочки ортогональных преобразований:

Объем арифметических действий алгоритма QL имеет порядок O(n³) на итерацию для матрицы общего вида, что является недопустимым. Однако для трехдиагональной матрицы этот объем составляет O(n) на итерацию, а для матрицы в форме Гессенберга O(n²), что делает подобные формы матриц высокоэффективными.

В этой секции мы рассматриваем случай действительной, симметричной, трехдиагональной матрицы. Таким образом, все собственные значения l_i действительны. Согласно теореме, если имеется l_j кратности p, то должно быть по крайней мере (p-1) нулей сверху и снизу от диагонали. Таким образом, матрицу можно будет разбить на подматрицы, которые можно будет диагонализировать раздельно, а усложнение диагональных блоков, имеющее место для матриц общего вида, здесь не важно.

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

Обобщая, представим, что уже найдены (r - 1) собственных значений A. Тогда матрицу можно уменьшить, обнулив первые (r - 1) столбцов и строк:
 

Значение k_s выбирается равным тому собственному значению ведущего минора размером (2 x 2), которое ближе к d_r. Можно показать, что скорость сходимости алгоритма при выборе этой стратегии обычно кубическая (и в худшем случае квадратичная при вырожденных собственных значениях). Эта исключительно быстрая сходимость делает алгоритм весьма привлекательным.

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

Как уже говорилось выше, QL-декомпозиция матрицы общего вида производится последовательностью преобразований Хаусхолдера. Однако для трехдиагональной матрицы значительно более эффективно применять плоские ротации Якоби P_pq. Элементы a₁₂, a₂₃, a₃₄, ..., a_n-1,n обнуляются серией ротаций P₁₂, P₂₃, P₃₄, ..., P_n-1,n. По симметрии, будут обнулены также поддиагональные элементы a₂₁, a₃₂, ..., a_n,n-1. Таким образом, каждая матрица трансформации Q_s будет представлять из себя произведение плоских ротаций:

Алгоритм основан на следующей лемме. Пусть A -- симметричная невырожденная матрица и B = Q^TAQ, где Q -- ортогональна, а B -- трехдиагональна с положительными внедиагональными элементами. Тогда Q и B полностью определяются на основании только последней строки Q^T. Доказательство. Пусть q_i^T обозначает i-ую строку матрицы Q^T, а q_i -- i-ый столбец матрицы Q. Соотношение BQ^T = Q^TA может быть записано как
 

n-ая строка этого матричного уравнения выглядит как

Для практического применения этой леммы предположим, что мы как-то нашли трехдиагональную матрицу A"_s+1, такую что

Далее, из оригинального алгоритма видно, что последняя строка Q_s^T совпадает с последней строкой P_n-1^(s). Но вспомним, что P_n-1^(s) служит для плоской ротации с целью обнуления элемента (n-1,n) матрицы A - k_s1. Простое вычисление показывает, что параметры c и s этой матрицы должны иметь вид:

Она должна быть приведена к трехдиагональной форме с помощью ортогональной матрицы, имеющей последнюю строку вида [0,0, ..., 0,1]; таким образом, последняя строка Q"_s^T останется равной строке P_n-1^(s). Это можно сделать последовательностью преобразований Хаусхолдера или Гивенса. Для матрицы указанного вида преобразование Гивенса эффективнее. Производим поворот в плоскости (n-2,n-1) для обнуления элемента (n-2,n) (по симметрии, будет также обнулен элемент (n, n-2)). После этого трехдиагональная форма сохраняется за исключением дополнительных элементов (n-3, n-1) и (n-1, n-3). Эти в свою очередь обнуляются поворотом в плоскости (n-3, n-2) и т.д. Таким образом, требуется последовательность (n - 2) трансформаций Гивенса. В результате будем иметь

Помещенная ниже программа tqli ("Tridiagonal QL Implicit"), алгоритмически основанная на опубликованных в [2, 3], на практике работает очень хорошо. Число итераций для первых собственных значений должно быть около 4 - 5, но помимо этого на этих итерациях также уменьшаются внедиагональные элементы правого нижнего угла. Таким образом, дальнейшие собственные значения выясняются значительно быстрее. Среднее число итераций на каждое собственное значение обычно составляет 1.3 - 1.6. Число действий на одной итерации имеет порядок O(n) с достаточно большим множителем, скажем ~20n. Таким образом, общее число действий по диагонализации оценивается как ~20n(1.5)n ~ 30n² операций. Если требуются собственные вектора, то число дополнительных операций будет значительно больше; их количество оценивается как 3n³.

Программа tqli

/* QL-алгоритм с неявными сдвигами для определения собственных значений (и собственных
   векторов) действительной, симметричной, трехдиагональной матрицы. Эта матрица может
   быть предварительно получена с помощью программы tred2. На входе d[1...n] содержит
   диагональ исходной матрицы, на выходе - собственные значения. На входе e[1...n]
   содержит поддиагональные элементы, начиная с e[2]. На выходе массив e разрушается.
   При необходимости поиска только собственных значений в программе следует
   закомментировать или удалить инструкции, необходимые только для поиска собственных
   векторов. Если требуются собственные вектора трехдиагональной матрицы, массив
   z[1...n][1...n] необходимо инициализировать на входе единичной матрицей. Если
   требуются собственные вектора матрицы, сведенной к трехдиагональному виду с помощью
   программы tred2, в массив z требуется загрузить соответствующий выход tred2. В
   обоих случаях на выходе массив z возвращает матрицу собственных векторов, расположенных
   по столбцам.
*/

/* максимальное число итераций */
#define MAXITER 30

void tqli(float *d, float *e, int n, float **z) {
  int m,l,iter,i,k;
  float s,r,p,g,f,dd,c,b;
  /* удобнее будет перенумеровать элементы e */
  for(i=2;i<=n;i++) e[i-1]=e[i];
  e[n]=0.;
  /* главный цикл идет по строкам матрицы */
  for(l=1;l<=n;l++) {
    /* обнуляем счетчик итераций для этой строки */
    iter=0;
    /* цикл проводится, пока минор 2х2 в левом верхнем углу начиная со строки l
       не станет диагональным */
    do {
      /* найти малый поддиагональный элемент, дабы расщепить матрицу */
      for(m=l;m<=n-1;m++) {
        dd=fabs(d[m])+fabs(d[m+1]);
        if((float)(fabs(e[m]+dd)==dd) break;
      }
      /* операции проводятся, если верхний левый угол 2х2 минора еще не диагональный */
      if(m!=l) {
        /* увеличить счетчик итераций и посмотреть, не слишком ли много. Функция
           nerror завершает программу с диагностикой ошибки. */
        if(++iter>=MAXITER) nerror("Too many iterations in tqli");
        /* сформировать сдвиг */
        g=(d[l+1]-d[l])/(2.*e[l]); r=hypot(1.,g);
        /* здесь d_m - k_s */
        if(g>=0.) g+=fabs(r);
        else g-=fabs(r);
        g=d[m]-d[l]+e[l]/g;
        /* инициализация s,c,p */
        s=c=1.; p=0.;
        /* плоская ротация оригинального QL алгоритма, сопровождаемая ротациями
           Гивенса для восстановления трехдиагональной формы */
        for(i=m-1;i>=l;i--) {
          f=s*e[i]; b=c*e[i];
          e[i+1]=r=hypot(f,g);
          /* что делать при малом или нулевом знаменателе */
          if(r==0.) {d[i+1]-=p; e[m]=0.; break;}
          /* основные действия на ротации */
          s=f/r; c=g/r; g=d[i+1]-p; r=(d[i]-g)*s+2.c*b; d[i+1]=g+(p=s*r); g=c*r-b;
          /* Содержимое следующего ниже цикла необходимо опустить, если
             не требуются значения собственных векторов */
          for(k=1;k<=n;k++) {
            f=z[k][i+1]; z[k][i+1]=s*z[k][i]+c*f; z[k][i]=c*z[k][i]-s*f;
          }
        }
        /* безусловный переход к новой итерации при нулевом знаменателе и недоведенной
           до конца последовательности ротаций */
        if(r==0. && i>=l) continue;
        /* новые значения на диагонали и "под ней" */
        d[l]-=p; e[l]=g; e[m]=0.;
      }
    } while(m!=l);
  }
}