Собственные вектора и значения матриц
Страница 8. Метод Хаусхолдера
Страница 8 из 10
 

Метод Хаусхолдера

Алгоритм Хаусхолдера приводит симметричную матрицу A размером (n x n) к трехдиагональной форме за (n - 2) ортогональных преобразований. Каждое преобразование уничтожает требуемую часть целой строки и целого соответствующего столбца. Базовой матрицей является матрица Хаусхолдера P, имеющая форму P = 1 - 2wwT,где w -- действительный вектор такой что |w|2 = 1. (В обозначениях, принятых здесь, матричное, или диадное, произведение векторов a и b записывается как abT, в то время как скалярное произведение этих векторов aTb). Матрица P является ортогональной, поскольку P2 = (1 - 2wwT)(1 - 2wwT) = 1 - 4wwT + 4w(wTw)wT = 1.Таким образом, P = P-1. Но PT = P, следовательно PT = P-1, что и доказывает ортогональность.

Перепишем P как

P = 1 - (uuT)/H,   где  H = |u|2/2,и теперь u может быть любым вектором. Пусть x -- вектор, составленный из первого столбца матрицы A. Выберем u = x - s|x|e1,где e1 -- единичный вектор [1, 0,...,0]T, s -- либо 1, либо -1, а выбор знака s будет сделан позже. Тогда Px = x - (u(x - s|x|e1)Tx)/H = x - (2u(|x|2 - s|x|x1))/(2|x|2 - 2s|x|x1) = x - u = s|x|e1.Это показывает, что матрица Хаусхолдера P действует на заданный вектор x, уничтожая все его элементы, кроме первого.

Для приведения симметричной матрицы A к трехдиагональной форме, вектор x, который будет образовывать первую матрицу Хаусхолдера, составим из нижних n - 1 элементов первого столбца. Тогда будут обнулены нижние n - 2 элемента:
 

P1A =  (
(
(
(
(
(
(		 1 0 0 ... 0 )
)
)
)
)
)
)		 x (
(
(
(
(
(
(		 a11 a12 a13 ... a1n )
)
)
)
)
)
)		  =  (
(
(
(
(
(
(		 a11 a12 a13 ... a1n )
)
)
)
)
)
)
0 (n-1)P1 a21 не важно k не важно                                    
0 a31 0                                          
... ... ...                                          
0 an1 0                                          

Здесь матрицы записаны в блочной форме, причем (n-1)P1 обозначает матрицу Хаусхолдера размерами ((n-1) x (n-1)). Величина k это просто плюс или минус модуль вектора [a21, ..., an1]T.

Полная ортогональная трансформация теперь будет выглядеть как:
 

A' = PAP =  (
(
(
(
(
(
(		 a11 k 0 ... 0 )
)
)
)
)
)
)
k не важно            
0              
...              
0              

Здесь был использован факт, что PT = P.

Для второй трансформации матрица Хаусхолдера составляется на основе нижних (n - 2) элементов второго столбца, и из нее конструируется матрица преобразования:
 

P2 =  (
(
(
(
(
(
(		 1 0 0 ... 0 )
)
)
)
)
)
)
0 1 0 ... 0      
0 0 (n-2)P2          
... ...            
0 0            

Единичный блок в левом верхнем углу обеспечивает сохранение трехдиагональности, полученной на предыдущей стадии, в то время как матрица Хаусхолдера (n-2)P2 размера (n-2) добавляет еще один столбец и строку к трехдиагональному результату. Ясно, что цепочка (n-2) таких преобразований приведет матрицу A к трехдиагональному виду.

Вместо того, чтобы явно проводить матричные умножения в выражениях PAP, мы вычислим вектор

p = (Au)/H.Тогда AP = A(1 - (uuT)/H) = A - puT,   A' = PAP = A - puT - upT + 2KuuT,где скаляр K определяется как: K = (uTp)/(2H).Если мы обозначим q = p - Ku,то будем иметь A' = A - quT - uqT.Последняя формула удобна для вычислений.

Программа редукции Хаусхолдера, основанная на алгоритме из [2] и помещенная ниже, реально начинает действие с n-ого столбца матрицы A, а не с первого, как было написано выше. В деталях алгоритм выглядит следующим образом. На стадии m (m = 1, 2, ..., n-2) вектор u имеет вид:

u = [ai1, ai2, ..., ai,i-2, ai,i-1 + s, 0, ..., 0].Здесь i = n - m + 1 = n, n-1, ..., 3.Величина s (что соответствует s|x| в прежних обозначениях) определяется как s2 = (ai1)2 + ... + (ai,i-1)2.Знак s выбирается так, чтобы совпадать со знаком ai,i-1 для уменьшения ошибки округления.

Переменные вычисляются в следующем порядке: s, u, H, p, K, q, A'. На каждой m-ой стадии матрица A становится трехдиагональной в последних m - 1 строках и столбцах.

Если уже найдены собственные вектора трехдиагональной матрицы (например, по программе из следующей секции), то собственные вектора A могут быть найдены умножением накопленного произведения трансформирующих матриц

Q = P1P2...Pn-2на матрицу этих векторов. Произведение трансформирующих матриц Q формируется рекурсией после нахождения всех Pi: Qn-2 = Pn-2, ..., Qj = PjQj+1, ..., Q = Q1.

На входе программы задается действительная симметричная матрица в массиве a[1...n][1...n]. На выходе этот массив содержит элементы матрицы Q. Вектор d[1...n] содержит множество диагональных элементов A', вектор e[1...n] содержит внедиагональные элементы A', начиная с e[2]. Заметим, что поскольку исходное содержимое массива a уничтожается, то при производстве серии последовательных вычислений с A он должен быть скопирован перед передачей в программу.

Промежуточные результаты не требуют дополнительной памяти. На стадии m вектора p и q не равны нулю только для элементов с индексами 1 ... i (вспомним, что i = n - m + 1), а вектор u -- только для элементов с индексами 1 ... i - 1. Элементы вектора e определяются в порядке n, n - 1, ..., так что вектор p можно держать на месте неопределенных элементов e. Вектор q может записываться поверх p, так как после определения q вектор p уже не требуется. Вектор u размещается в i-ой строке массива a, а u/H -- в i-ом столбце. После проведения редукции, матрицы Qj вычисляются с использованием записанных в массив a векторов u и u/H. Поскольку матрица Qj является единичной для последних n - j + 1 строк и столбцов, то требуется вычислить ее элементы лишь до строки и столбца n - j. При этом можно перезаписать u и u/H в соответствующих строке и столбце, поскольку они более не нужны.

Расположенная ниже программа tred2 содержит еще одну тонкость. Если величина s2 равна нулю или "мала" на какой-либо из стадий, то соответствующее преобразование может быть пропущено. Простой критерий типа

s2 < (наименьшее представимое число)/(машинная точность)в большинстве случаев является слишком сильным. В действительности берется более аккуратный критерий. Определим число e = Sk=1i-1(|aik|).При e от нуля до машинной точности преобразование пропускается. В противном случае мы переопределяем aik --> aik/eи проводим преобразование с масштабированными величинами. (Преобразование Хаусхолдера зависит только от отношения элементов).

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

Программа tred2 предназначена для совместного использования с программой tqli, помещенной в следующей секции. tqli находит собственные вектора и значения симметричной трехдиагональной матрицы. Комбинация tred2 и tqli на сегодня является наиболее эффективным алгоритмом поиска всех собственных чисел (и, при необходимости, собственных векторов) действительной симметричной матрицы.

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

Программа tred2

/* Редукция Хаусхолдера действительной симметричной матрицы a[1...n][1...n].
   На выходе a заменяется ортогональной матрицей трансформации q.
   d[1...n] возвращает диагональ трехдиагональной матрицы.
   e[1...n] возвращает внедиагональные элементы, причем e[1]=0.
   Некоторые инструкции программы могут быть опущены (это указано в комментариях),
   если требуется отыскать только собственные значения, а не вектора. В
   этом случае на выходе массив a будет содержать мусор.
*/
#include <math.h>

void tred2(float **a, int n, float *d, float *e) {
  int l,k,j,i;
  float scale,hh,h,g,f;
  /* Проход по стадиям процесса редукции */
  for(i=n;i>=2;i--) {
    l=i-1; h=scale=0.;
    /* сложный процесс везде, кроме последней стадии */
    if(l>1) {
      /* вычислить шкалу */
      for(k=1;k<=l;k++) scale += fabs(a[i][k]);
      /* малая величина шкалы -> пропустить преобразование */
      if(scale==0.) e[i]=a[i][l];
      else {
       /* отмасштабировать строку и вычислить s2 в h */
        for(k=1;k<=l;k++) {
          a[i][k]/=scale; h += a[i][k]*a[i][k];
        }
       /* вычислить вектор u */
        f=a[i][l];
        g=(f>=0.?-sqrt(h):sqrt(h));
        e[i]=scale*g; h -= f*g;
       /* записать u на место i-го ряда a */
        a[i][l]=f-g;
       /* вычисление u/h, Au, p, K */
        f=0.;
        for(j=1;j<=l;j++) {
       /* следующая инструкция не нужна, если не требуются вектора,
             она содержит загрузку u/h в столбец a */
          a[j][i]=a[i][j]/h;
       /* сформировать элемент Au (в g) */
          g=0.;
          for(k=1;k<=j;k++) g += a[j][k]*a[i][k];
          for(k=j+1;k<=l;k++) g += a[k][j]*a[i][k];
       /* загрузить элемент p во временно неиспользуемую область e */
          e[j]=g/h;
       /* подготовка к формированию K */
          f += e[j]*a[i][j];
        }
       /* Сформировать K */
        hh=f/(h+h);
        for(j=1;j<=l;j++) {
       /* Сформировать q и поместить на место p (в e) */
          f=a[i][j]; e[j]=g=e[j]-hh*f;
       /* Трансформировать матрицу a */
          for(k=1;k<=j;k++) a[j][k] -= (f*e[k]+g*a[i][k]);
        }
      }
    }
    else e[i]=a[i][l];
    d[i]=h;
  }
  /* если не нужны собственные вектора, опустите следующую инструкцию */
  d[1]=0.;
  /* эту опускать не надо */
  e[1]=0.;
  /* Все содержание цикла, кроме одной инструкции, можно опустить, если не
     требуются собственные вектора */
  for(i=1;i<=n;i++) {
    l=i-1;
    /* этот блок будет пропущен при i=1 */
    if(d[i]!=0.) {
      for(j=1;j<=l;j++) {
        g=0.;
       /* формируем PQ, используя u и u/H */
        for(k=1;k<=l;k++) g += a[i][k]*a[k][j];
        for(k=1;k<=l;k++) a[k][j] -= g*a[k][i];
      }
    }
    /* эта инструкция остается */
    d[i]=a[i][i];
    /* ряд и колонка матрицы a преобразуются к единичной, для след. итерации */
    a[i][i]=0.;
    for(j=1;j<=l;j++) a[j][i]=a[i][j]=0.;
  }
}
 
« Предыдущая статья   Следующая статья »