La factorización QR se utiliza para descomponer una matriz (que puede representar un sistema de ecuaciones a resolver) en dos matrices, con las cuales es más fácil de operar y hallar soluciones al sistema de ecuaciones.
Más información en la WIkipedia.
Los siguientes códigos están basados en los algoritmos mostrados en el libro "Numerical Linear Algebra" de Trefethen y Bau.
Algoritmo de Gram-Schmidth Modificado (devuelve la factorización QR reducida):
El código en Matlab es:
function [Q, R] = gs(A)
[m, n] = size(A);
R = zeros(n, n);
V = A;
Q=zeros(m, n);
for i =1:n
R(i,i)=
norm(V(:,i));
Q(:,i)=
V(:,i)/R(i,i);
for j=i+1:n
R(i,j)=
(Q(:,i)')*V(:,j);
V(:,j)=V(:,j)
- R(i,j)*Q(:,i);
end
end
"A" es una matriz de rango completo, la función devuelve dos matrices: Q (la matriz ortonormal) y R (la matriz triangular).
Problemas con este algoritmo: Sólo da la factorización QR reducida, y sólo funciona para matrices de
rango completo (es decir: ninguna columna es múltiplo de, y/o resultado de, sumar otras columnas de la mariz entre sí).
Modifiqué este algoritmo para que funcionara también con matrices de "rango deficiente" (lo opuesto a "rango completo"), y que devolviera la factorización completa; esto último se hace completando las columnas de la matriz Q, las cuales son ortogonales entre sí. En palabras del curso de Álgebra Lineal: "completando la base ortonormal del espacio vectorial generado por las columnas de Q". Dicho de otra forma: aplicando el
Teorema de la Base Ortonormal Incompleta.
El código en Matlab es:
function [Q, R] = gsC(A)
clc;
[m, n] = size(A);
R = zeros(m, n);
V = A;
Q=zeros(m, n);
k=1; %
inicializo eta variable, sólo la uso para matrices de rango deficiente
for i =1:n
R(i,i)=
norm(V(:,i));
if abs(R(i,i))
> 1.0e-10
Q(:,i)=
V(:,i)/R(i,i);
else % salió cero! La matriz A es de
rango deficiente
r= zeros(1, m); % este va
a ser el vector que reemplazará a la fila problemática de A
Q(:,i)= r;
Q(k,
i)=1;
k=k+1;
Q2 = Q; % Q2 tendrá los valores de
Q, en Q irá la nueva base ortonormal
% ortogonalizo otra vez aplicando
Gram-Schmidt a la matriz Q para
% hallar un vector con el qué
reemplazar a la fila problemática y
% seguir teniendo una base ortonormal
for j = 1 : i-1
r(j)
= (Q2(:,j)')*Q2(:,i);
Q(:,i)=Q(:,i)
- r(j)*Q2(:,j);
Q(:,i)=Q(:,i)/norm(Q(:,i));
end
end
for j=i+1:n
R(i,j)=
(Q(:,i)')*V(:,j);
V(:,j)=V(:,j)
- R(i,j)*Q(:,i);
end
end
% hasta
aquí Q y R son la descomposición reducida
% acá
hallo la completa, Gram Schmidt otra vez
if m > n
I = eye(m,
m - n);
Q=[Q I];
V=Q;
r= zeros(1, m);
k=0;
for i=n+1:m
for j = 1 : n + k
r(j)
= (V(:,j)')*V(:,i);
Q(:,i)=Q(:,i)
- r(j)*V(:,j);
Q(:,i)=Q(:,i)/norm(Q(:,i));
end
k=k+1;
end
end
"A" es una matriz cualquiera, la función devuelve dos matrices: Q (la matriz ortonormal) y R (la matriz triangular).
El código en Matlab es:
function [Q, R] =
houseHolder(A)
clc;
[m n] = size(A);
R= A;
Q = eye(m,m);
n1=min(m,n);
for i= 1:n1
x = R(i:m,
i);
s=0;
if x(1)==0
s = 1;
else
s=
sign(x(1));
end
e = zeros(m
- i + 1, 1);
e(1)=1;
u = s * norm(x) * e + x;
u = u / norm(u);
R(i:m,
i:n) = R(i:m, i:n) - 2*(u*u')*R(i:m,
i:n);
Q(:,i:m) =
Q(:,i:m) - Q(:,i:m)*(2*u*u');
end
Igual que en los casos anteriores, "A" es una matriz cualquiera, la función devuelve dos matrices: Q (la matriz ortonormal) y R (la matriz triangular).
Los tres algoritmos han sido probados en Matlab 2008a, así que deben correr en cualquier versión posterior, y no tan posterior :)
