Домівка > Чисельні методи > Порівняльний тест стабільності процесу Грама — Шмідта

Порівняльний тест стабільності процесу Грама — Шмідта

Сьогодні ми розглянемо стабільність процесу Грама-Шмідта. Для цього ми реалізуємо класичну, clgs(), і стабільну, stgs(), версії алгоритму, а як еталон ми використаємо вбудовану функцію MATLAB – qr(). Ми згенеруємо 80×80 матрицю з “випадковою” Q і з R чиї діагональні елементи спадають експоненційно. Для створення вхідної матриці A ми скористаємось cингулярним розкладом матриці. Більше про те чим відрізняєтьс класичний і стабільний варіанти алгоритму можна прочитати у статті на Вікіпедії.

function [Q,R] = clgs(A)
  Q = zeros(size(A));
  R = zeros(size(A));
  
  a1 = A(:,1);
  R(1,1) = sqrt(dot(a1,a1));
  Q(:,1) = a1 / R(1,1);
  
  for i = 2:size(A,2)
    v = A(:,i);
    p = zeros(size(A,2),1);
    for j = 1:i-1
      u = Q(:,j); % u - already unit length
      R(j,i) = dot(v,u);
      p = p + R(j,i) * u;
    end
    v = v - p;
    R(i,i) = sqrt(dot(v,v));
    Q(:,i) = v / R(i,i);
  end
end
function [Q,R] = stgs(A)
  Q = zeros(size(A));
  R = zeros(size(A));
  
  a1 = A(:,1);
  R(1,1) = sqrt(dot(a1,a1));
  Q(:,1) = a1 / R(1,1);
  
  for i = 2:size(A,2)
    v = A(:,i);
    for j = 1:i-1
      u = Q(:,j); % u - already unit length
      R(j,i) = dot(v,u);
      p = R(j,i) * u;
      v = v - p;
    end
    R(i,i) = sqrt(dot(v,v));
    Q(:,i) = v / R(i,i);
  end
end

Як бачите відмінності між clgs і stgs дуже маленькі. В stgs ми оновлюємо v на кожній ітерації внутрішнього циклу і використовуємо його для обчислення наступної проекції.

[U,X] = qr(randn(80)); % randn() використовує нормальний розподіл.
[V,X] = qr(randn(80));
S = diag(2 .^ (-1:-1:-80)); % .^ означає поелементне піднесення до степеня.
A = U*S*V;

[Qc,Rc] = clgs(A); % Тут нас цікавить тільки R, але
[Qs,Rs] = stgs(A); % також можна перевірити norm(Q’*Q-eye(80)).
[Q,R] = qr(A);
x = diag(Rc); % Так ми збираємо діагональні елементи у вектор.
y = diag(Rs);
z = abs(diag(R)); % qr() не гарантує додатну діагональ.

abscissae=1:1:80;
plot(abscissae,log2(x),'+',abscissae,log2(y),'o',abscissae,log2(z),'*');
legend('clgs', 'stgs', 'qr');

Ми обрали U і V випадковим чином, а S має діагональні елементи, що спадають експоненційно. Так ми отримаємо діагональні елементи R приблизно того ж порядку. Щонайменше деякі з них загублятся через похибку заокруглення.


Зауважте, що похибка з’їла clgs близько 2−28, і це задовго до того як машинне епсілон мало стати проблемою. (\epsilon_{machine} = 2^{-56} = (2^{-28})^2.) Тоді як, обидва stgs і qr залишаються стабільними аж до границі машинного епсілон. Чи можемо ми їх розрізнити іншим способом?

Advertisements
  1. Коментарів ще немає.
  1. No trackbacks yet.

Залишити відповідь

Заповніть поля нижче або авторизуйтесь клікнувши по іконці

Лого WordPress.com

Ви коментуєте, використовуючи свій обліковий запис WordPress.com. Log Out / Змінити )

Twitter picture

Ви коментуєте, використовуючи свій обліковий запис Twitter. Log Out / Змінити )

Facebook photo

Ви коментуєте, використовуючи свій обліковий запис Facebook. Log Out / Змінити )

Google+ photo

Ви коментуєте, використовуючи свій обліковий запис Google+. Log Out / Змінити )

З’єднання з %s

%d блогерам подобається це: