Acerca de:

Este blog contiene los códigos, ejemplos y bases de datos que he usado cuando aprendía acerca de algún tema específico. En lugar de borrarlos (una vez dominado ya el tema), he decidido publicarlos :)

domingo, 24 de agosto de 2014

Algoritmos para calcular la factorización QR de una matriz

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). 

Otra forma de calcular la factorización QR (y la que utiliza internamente Matlab) es usando la factorización de Householder:


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 :)