Сведение симметричной матрицы к трехдиагональной форме: редукции Гивенса и Хаусхолдера

Как уже было сказано, оптимальная стратегия для поиска собственных значений и собственных векторов заключается в первоначальном сведении матрицы к возможно более простой форме с дальнейшем применением итерационной процедуры. Для симметричных матриц наиболее предпочтительной простейшей формой является трехдиагональная. Редукция Гивенса является модификацией метода Якоби. Вместо того, чтобы стараться все время приводить матрицу к диагональному виду, остановка в этом методе происходит, когда матрица трехдиагональна. Это позволяет процедуре завершиться за конечное число шагов в отличие от метода Якоби, требующего последовательные итерации для сходимости.

Метод Гивенса

При использовании метода Гивенса угол ротации выбирается таким, чтобы обнулить элемент, находящийся не на одном из четырех углов, т.е. не app, apq или aqq. Конкретно, матрица P23 используется для обнуления a31 (и по симметрии a13). Затем P24 используется для обнуления a41. В целом, используется последовательность
P23, P24, ..., P2n; P34, ..., P3n; ...; Pn-1,n
для обнуления ak,j-1 с помощью Pjk. Метод работает, поскольку такие элементы, как a'rp и a'rq при r!=p и r!=q суть линейные комбинации старых значений arp и arq. Таким образом, если arp и arq уже нулевые, то они таковыми останутся и после ротации. Очевидно, требуется порядка n2/2 ротаций, а число умножений в алгоритме будет порядка 4n3/3, не считая тех, которые требуются для преобразования произведения трансформирующих матриц, необходимого для вычисления собственных векторов.

Метод Хаусхолдера, который будет описан ниже, является таким же устойчивым, как и редукция Гивенса, но требует в 2 раза меньше операций, таким образом, метод Гивенса обычно не используется.  Однако недавние исследования [6] показали, что редукция Гивенса может быть переформулирована так, чтобы повысить ее эффективность примерно в 2 раза, а также избежать вычисления квадратных корней. В результате этого "быстрый метод Гивенса" становится конкурентоспособным с методом Хаусхолдера. С другой стороны, быстрый метод должен отслеживать возможности переполнения, и переменные должны быть периодически перемасштабируемы при его использовании. Это обстоятельство говорит о том, что метод Хаусхолдера все же предпочтительнее.

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

Алгоритм Хаусхолдера приводит симметричную матрицу 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.;
  }
}

К оглавлению
Hosted by uCoz