|
Страница 6 из 10
Преобразования Якоби симметричной матрицы Метод Якоби состоит из цепочки ортогональных преобразований подобия матрицы. Каждое преобразование (ротация Якоби) -- это плоский поворот с целью обнуления одного из внедиагональных элементов матрицы. Последовательные преобразования не сохраняют уже установленные нулевые элементы, но вместе с тем внедиагональные элементы становятся меньше и меньше до тех пор, пока матрица не станет диагональной с точностью до машинного нуля. Накопление в процессе преобразований произведения трансформационных матриц дает матрицу собственных векторов, в то время как диагональные элементы являются собственными значениями. Метод Якоби всегда является робустным для действительных симметричных матриц. Для матриц размера, скажем больше 10, этот алгоритм более медленный, чем метод QR (см. ниже). Однако алгоритм Якоби значительно проще других более эффективных методов, таким образом он рекомендуется в случаях, когда время выполнения не является критическим. Базовая матрица ротации Якоби имеет вид:
| Ppq | = | (
(
(
(
(
(
(
(
( | 1 | | | | | | | )
)
)
)
)
)
)
)
) | | | ... | | | | | | | | | | | | | c | ... | s | | | | | | | | | | ... | 1 | ... | | | | | | | | | | -s | ... | c | | | | | | | | | | | | | ... | | | | | | | | | | | | | 1 | | | | | Все диагональные элементы матрицы ротации Якоби равны 1, за исключением двух элементов в строках (и столбцах) p и q. Все внедиагональные элементы равны нулю, за исключением двух элементов, равных s и (-s). Числа c и s являются косинусом и синусом угла поворота f, так что c2+s2=1. Ротация Якоби преобразует матрицу A согласно формуле: A' = PpqTAPpq.Произведение PpqTA изменяет только строки p и q матрицы A, в то время как APpq меняет только столбцы p и q. Отметим, что нижние индексы p и q в обозначении Ppq не означают элемент этой матрицы, а скорее индексируют тип ротации, т.е. плоскость, в которой она происходит. Таким образом, измененные элементы в матрице A' расположены только в строках и столбцах p и q:
| A' | = | (
(
(
(
(
(
(
(
( | | ... | a'1p | ... | a'1q | ... | | )
)
)
)
)
)
)
)
) | | ... | | ... | | ... | | ... | | | | | | a'p1 | ... | a'pp | ... | a'pq | ... | a'pn | | | | | | ... | | ... | | ... | | ... | | | | | | a'q1 | ... | a'qp | ... | a'qq | ... | a'qn | | | | | | | | ... | | ... | | ... | | | | | | | ... | a'np | ... | a'nq | ... | | | | | | Для измененных элементов матрицы (с учетом приведенного выше и ее симметрии) получаются явные формулы: a'rp = carp - sarq (r!=p, r!=q)
a'rq = carq + sarp (r!=p, r!=q)
a'pp = c2app + s2aqq - 2csapq
a'qq = s2app + c2aqq + 2csapq
a'pq = (c2 - s2)apq + cs(app- aqq)Идея метода Якоби -- в том, чтобы постараться обнулить внедиагональные элементы серией ротаций. Для того, чтобы обнулилось a'pq, угол ротации, согласно формулам, приведенным выше, должен быть следующим: q = cot(2f) = (c2 - s2)/(2cs) = (aqq- app)/(2apq).Полагая t=s/c, для определения t прийдем к следующему уравнению: t2 + 2tq - 1 = 0.Меньший по модулю корень этого уравнения соответствует вращению на угол, меньший p/4; выбор этого угла на каждой из стадий дает наиболее устойчивый алгоритм. Используя формулу для решения квадратного уравнения с дискриминантом в знаменателе, меньший корень определяется как: t = sign(q)/(|q|+sqrt(q2+1)).Если величина q такова, что q2 даст переполнение, полагаем t = 1/(2q). По значению t определяются c и s: c = 1/sqrt(t2+1), s = tc.При численной модификации элементов матрицы мы стремимся к уменьшению ошибки округления. Идея состоит в том, чтобы переписать их как прежнее значение плюс малая корректирующая добавка. В этом случае коэффициенты матрицы после преобразования будут выглядеть как: a'pp = app - tapq
a'rp = arp - s(arq + garp)
a'rq = arq + s(arp - garq) ,где t (тангенс половины угла поворота) определяется, как t = s/(1+c). Сходимость метода Якоби можно оценить, рассматривая сумму квадратов внедиагональных элементов S; уравнения преобразований дают для этой суммы после трансформации следующее соотношение: S = Sr!=s(|ars|2); S' = S - 2|apq|2.Поскольку преобразования ортогональные, сумма диагональных элементов при этом возрастает соответственно на |apq|2. Сумма же квадратов внедиагональных элементов монотонно убывает. Поскольку она ограничена снизу нулем и поскольку мы можем каждый раз выбирать какой хотим элемент apq для обнуления, сумма будет стремиться к нулю. После цепочки преобразований получается матрица D, диагональная с точностью до машинного нуля. Ее диагональные элементы образуют собственные значения первоначальной матрицы A, поскольку выполняется D = VTAV, где V = P1P2P3...,а Pi -- матрицы ротации Якоби. Столбцы матрицы V являются собственными векторами. Они вычисляются рекуррентно, на каждой стадии по формулам: V' = VPi ,где первоначально V -- единичная матрица. Для преобразования элементов V используются формулы: v'rs = vrs (s!=p, s!=q)
v'rp = cvrp - svrq
v'rq = svrp + cvrqДля уменьшения ошибки округления вышеприведенные формулы следует переписать в терминах g (как выше для элементов матрицы A). Единственным остающимся вопросом теперь является стратегия, по которой нужно перебирать уничтожаемые элементы. Оригинальный алгоритм Якоби от 1846 года на каждой стадии искал в верхней треугольной области наибольший по модулю элемент и обнулял его. Это весьма рациональная стратегия для ручного счета, но она неприемлема для компьютера, поскольку делает число операций на каждой стадии ротации Якоби порядка N2 вместо N. Для наших целей лучшей стратегией является циклический метод Якоби, где элементы удаляются в строгом порядке. Например, можно просто проходить по строкам: P12 , P13 , ... P1n ; P23 ,...,P2n ;... Можно показать, что скорость сходимости будет в общем случае квадратичной, как для оригинального, так и для циклического метода Якоби. Назовем одно такое множество последовательных n(n-1)/2 преобразований Якоби проходом. Программа, помещенная ниже, основана на алгоритмах из [2, 3]; она содержит еще две тонкости: - На первых трех проходах мы проводим ротацию относительно индекса (pq) только тогда, когда |apq|>e для некоторого порогового значения e = S0/(5n2), где S0 -- сумма модулей внедиагональных элементов верхней треугольной матрицы.
- После четырех проходов в том случае, когда |apq|<<|app| и |apq|<<|aqq|, полагаем |apq|=0 и пропускаем ротацию. Критерием сравнения является превышение диагональных элементов в 10D+2 раз, где D -- число значащих десятичных цифр.
В приведенной ниже программе симметричная матрица a размером (n x n) на входе загружена в массив a[1...n][1...n]. На выходе наддиагональные элементы a разрушаются, но диагональные и поддиагональные остаются на месте с тем, чтобы сохранить информацию об исходной симметричной матрице. Вектор d[1...n] возвращает собственные значения a. Во время вычисления он содержит текущую диагональ. Матрица v[1...n][1...n] возвращает набор нормализованных собственных векторов, относящихся к d[k] для ее k-го столбца. Параметр nrot -- число ротаций Якоби, необходимых для достижения сходимости. Типичные матрицы требуют от 6 до 10 проходов для достижения сходимости, или от 3n2 до 5n2 ротаций. Каждая ротация требует порядка 4n операций умножения и сложения, так что общие затраты имеют порядок от 12n3 до 20n3 операций. При вычислении собственных векторов наряду с собственными значениями, число операций на ротации возрастает с 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 или других программ поиска собственных векторов и значений (например, расположенных ниже). Метод сортировки -- прямая вставка -- имеет порядок N2 (значительно больше, чем Nlog(N) для эффективных методов, но по сравнению с порядком N3 для процедуры 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;
}
}
} |