Собственные вектора и значения матриц
Страница 10. QL-алгоритм с неявными сдвигами
Страница 10 из 10
 

QL-алгоритм с неявными сдвигами

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

Алгоритм основан на следующей лемме. Пусть A -- симметричная невырожденная матрица и B = QTAQ, где Q -- ортогональна, а B -- трехдиагональна с положительными внедиагональными элементами. Тогда Q и B полностью определяются на основании только последней строки QT. Доказательство. Пусть qiT обозначает i-ую строку матрицы QT, а qi -- i-ый столбец матрицы Q. Соотношение BQT = QTA может быть записано как
 

(
(
(
(
(
(
(
(		 b1 g1           )
)
)
)
)
)
)
)		  * (
(
(
(
(
(
(
(		 q1T )
)
)
)
)
)
)
)		  =  (
(
(
(
(
(
(
(		 q1T )
)
)
)
)
)
)
)		 A
a2 b2 g2         q2T q2T                  
      ...       ... ...                  
        an-1 bn-1 gn-1 qn-1T qn-1T                  
          an bn qnT qnT                  

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

anqn-1T +  bnqnT = qnTA.Поскольку матрица Q ортогональна, qnTqm = dnm.Если теперь умножим обе части уравнения на qn, то найдем: bn = qnTAqn,которая однозначно известна при заданном qn. Кроме того, из уравнения следует: anqn-1T = zn-1T,  где обозначено zn-1T = qnTA - bnqnT,последняя величина также известна. Таким образом, an2 = zn-1Tzn-1,  или an = |zn-1|, а также qn-1T =  zn-1T/an.(Поскольку согласно гипотезе an ненулевое). Аналогично можно по индукции доказать, что если нам известны qn, qn-1, ..., qn-j, а также все a, b и g вплоть до уровня (n - j), то мы можем определить все значения на уровне (n - (j+1)).

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

A"s+1 = Q"sTA"sQ"s,где Q"sT ортогональна и ее последняя строка совпадает с последней строкой матрицы QsT из оригинального QL-алгоритма. Тогда будем иметь Q"s = Qs и A"s+1 = As+1.

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

c = (dn - ks)/sqrt(en2 + (dn - ks)2),  s = (- en-1)/sqrt(en2 + (dn - ks)2).Матрица Pn-1(s)AsPn-1(s)T является трехдиагональной с двумя дополнительными элементами:
 
(
(
(
(
(
(
(
(		 ...           )
)
)
)
)
)
)
)
  ...            
  x x x        
    x x x x    
      x x x    
      x x x    

Она должна быть приведена к трехдиагональной форме с помощью ортогональной матрицы, имеющей последнюю строку вида [0,0, ..., 0,1]; таким образом, последняя строка Q"sT останется равной строке Pn-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) трансформаций Гивенса. В результате будем иметь

QsT = Q"sT = P"1(s)P"2(s)...P"n-2(s)Pn-1(s),где P" -- матрицы трансформаций Гивенса, а Pn-1 -- плоская ротация из оригинального алгоритма. По значению QsT находим значение A на новой итерации. Заметим, что сдвиг ks производится неявно при вычислении параметров c и s.

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

Программа tqli

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

/* 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);
  }
}

Литература по проблеме собственных векторов и значений

[1] Stoer, J., and Burlisch, R. 1980. Introduction to Numerical Analysis (New York: Springer Verlag, Chapter 6).
[2] Wilkinson, J.H., and Reinisch, C. 1971. Linear Algebra, vol. ii of Handbook for Automatic Computation (New York: Springer Verlag).
[3] Smith, B.T., et al. 1976. Matrix Eigensystem Routines -- EISPACK Guide, 2nd ed., vol. 6 of Lecture Notes in Computer Science (New York, Springer Verlag).
[4] IMSL Math/Library Users Manual (IMSL Inc., 2500 CytyWest Boulevard, Houston, TX 77042).
[5] NAG Fortran Library (Numerical Algorithms Group, 256 Banbury Road, Oxford OX27DE, U.K.).
[6] Golub, G.H., and Van Loan, C.F. 1989. Matrix Computations, 2nd ed. (Baltimore: Johns Hopkins University Press), 5.1, 7.7, 8.4.
[7] Wilkinson, J.H. 1965. The Algebrain Eigenvalue Problem (New York: Oxford University Press).
[8] Acton, F.S. 1970. Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), Chapter 13.
[9] Horn, R.A., and Johnson, C.R. 1985. Matrix Analysis (Cambridge: Cambridge University Press).

« Предыдущая страница - Следующая страница
 
« Предыдущая статья   Следующая статья »