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