12.VL: QR-Zerlegung, Gram-Schmidt Verfahren

1. Test zur Genauigkeit

Called Functions

gramschmidt
modgramschmidt

d=1.d-8;
A = [1,1,1;d,0,0;0,d,0;0,0,d];
[Q,R] = gramschmidt(A);
Q'*Q

ans =

    1.0000   -0.0000   -0.0000
   -0.0000    1.0000    0.5000
   -0.0000    0.5000    1.0000

[Q,R] = modgramschmidt(A);
Q'*Q

ans =

    1.0000   -0.0000   -0.0000
   -0.0000    1.0000    0.0000
   -0.0000    0.0000    1.0000

2. Test zur Genauigkeit

dims = [1:160];
err = 0* dims;
errmod = 0*dims;
errref = 0*dims;

for dim = dims
  A = randn(dim);
  [Q,R] = gramschmidt(A);
  err(dim) = max(max(abs(Q*Q'-eye(dim))));
  [Q,R] = modgramschmidt(A);
  errmod(dim) = max(max(abs(Q*Q'-eye(dim))));
  [Q,R] = qr(A);
  errref(dim) = max(max(abs(Q*Q'-eye(dim))));
end

semilogy(dims,err,'x',dims,errmod,'o',dims,errref,'r.')
 shg

% Gram-Schmidt Orthogonalisierung
function [Q,R] = gramschmidt(A)
% Eingabe: Matrix A (n x m)
% Ausgabe: Matrix Q (n x m) mit orthonormalen Spalten
%          Matrix R (m x m) obere Dreiecksmatrix

[n,m] = size(A);

Q = zeros(n,m);
R = zeros(m,m);
for j = 1:m
   Q(:,j) = A(:,j);
   % Orthogonalisiere A(:,j) gegen Q(:,1),...,Q(:,j-1)
   for i = 1:j-1
      R(i,j) = Q(:,i)'*A(:,j);
      Q(:,j) = Q(:,j) - R(i,j)*Q(:,i);
   end
   % Normiere Q(:,j)
   R(j,j) = norm(Q(:,j));
   Q(:,j) = Q(:,j)/R(j,j);
end

% Modifizierte Gram-Schmidt Orthogonalisierung
function [Q,R] = modgramschmidt(A)
% Eingabe: Matrix A (n x m)
% Ausgabe: Matrix Q (n x m) mit orthonormalen Spalten
%          Matrix R (m x m) obere Dreiecksmatrix

[n,m] = size(A);

Q = zeros(n,m);
R = zeros(m,m);
for j = 1:m
   Q(:,j) = A(:,j);
   % Orthogonalisiere Q(:,j) gegen Q(:,1),...,Q(:,j-1)
   for i = 1:j-1
      % Berechne die Koeffizienten R(i,j) mit aktuellem Q(:,j) statt A(:,j)
      R(i,j) = Q(:,i)'*Q(:,j);
      Q(:,j) = Q(:,j) - R(i,j)*Q(:,i);
   end
   % Normiere Q(:,j)
   R(j,j) = norm(Q(:,j));
   Q(:,j) = Q(:,j)/R(j,j);
end