Factorización QR

Recurriendo al método de ortogonalización de Gram-Schmidt, con las columnas de A como los vectores a procesar.

${\displaystyle A=(\mathbf {a} _{1}|\cdots |\mathbf {a} _{n})}$. Entonces

${\displaystyle \mathbf {u} _{1}=\mathbf {a} _{1},\qquad \mathbf {e} _{1}={\mathbf {u} _{1} \over \|\mathbf {u} _{1}\|}}$

${\displaystyle \mathbf {u} _{2}=\mathbf {a} _{2}-\mathrm {proj} _{\mathbf {e} _{1}}\,\mathbf {a} _{2},\qquad \mathbf {e} _{2}={\mathbf {u} _{2} \over \|\mathbf {u} _{2}\|}}$

${\displaystyle \mathbf {u} _{3}=\mathbf {a} _{3}-\mathrm {proj} _{\mathbf {e} _{1}}\,\mathbf {a} _{3}-\mathrm {proj} _{\mathbf {e} _{2}}\,\mathbf {a} _{3},\qquad \mathbf {e} _{3}={\mathbf {u} _{3} \over \|\mathbf {u} _{3}\|}}$

${\displaystyle \vdots }$

${\displaystyle \mathbf {u} _{k}=\mathbf {a} _{k}-\sum _{j=1}^{k-1}\mathrm {proj} _{\mathbf {e} _{j}}\,\mathbf {a} _{k},\qquad \mathbf {e} _{k}={\mathbf {u} _{k} \over \|\mathbf {u} _{k}\|}}$

Naturalmente, utilizamos los a_is de A para obtener:

${\displaystyle \mathbf {a} _{1}=\mathbf {e} _{1}\|\mathbf {u} _{1}\|}$

${\displaystyle \mathbf {a} _{2}=\mathrm {proj} _{\mathbf {e} _{1}}\,\mathbf {a} _{2}+\mathbf {e} _{2}\|\mathbf {u} _{2}\|}$

${\displaystyle \mathbf {a} _{3}=\mathrm {proj} _{\mathbf {e} _{1}}\,\mathbf {a} _{3}+\mathrm {proj} _{\mathbf {e} _{2}}\,\mathbf {a} _{3}+\mathbf {e} _{3}\|\mathbf {u} _{3}\|}$

${\displaystyle \vdots }$

${\displaystyle \mathbf {a} _{k}=\sum _{j=1}^{k-1}\mathrm {proj} _{\mathbf {e} _{j}}\,\mathbf {a} _{k}+\mathbf {e} _{k}\|\mathbf {u} _{k}\|}$

Ahora estas ecuaciones pueden ser escritas en forma matricial de esta manera:

${\displaystyle \left(\mathbf {e} _{1}\left|\ldots \right|\mathbf {e} _{n}\right){\begin{pmatrix}\|\mathbf {u} _{1}\|&\langle \mathbf {e} _{1},\mathbf {a} _{2}\rangle &\langle \mathbf {e} _{1},\mathbf {a} _{3}\rangle &\ldots \\0&\|\mathbf {u} _{2}\|&\langle \mathbf {e} _{2},\mathbf {a} _{3}\rangle &\ldots \\0&0&\|\mathbf {u} _{3}\|&\ldots \\\vdots &\vdots &\vdots &\ddots \end{pmatrix}}}$ :::::::::

El producto de cada fila con cada columna de las matrices de arriba, nos da la respectiva columna de A con la que comenzamos y, por tanto, dada la matriz A, la hemos factorizado en una matriz ortogonal Q (la matriz de e_ks), aplicando el proceso de Gram-Schmidt, y la matriz resultante triangular superior es R.

Alternativamente, la matriz ${\displaystyle {\begin{matrix}R\end{matrix}}}$ puede calcularse de la siguiente manera:

Recordemos que: ${\displaystyle {\begin{matrix}Q\end{matrix}}=\left(\mathbf {e} _{1}\left|\ldots \right|\mathbf {e} _{n}\right).}$ Entonces, tenemos

${\displaystyle {\begin{matrix}R=Q^{T}A=\end{matrix}}{\begin{pmatrix}\langle \mathbf {e} _{1},\mathbf {a} _{1}\rangle &\langle \mathbf {e} _{1},\mathbf {a} _{2}\rangle &\langle \mathbf {e} _{1},\mathbf {a} _{3}\rangle &\ldots \\0&\langle \mathbf {e} _{2},\mathbf {a} _{2}\rangle &\langle \mathbf {e} _{2},\mathbf {a} _{3}\rangle &\ldots \\0&0&\langle \mathbf {e} _{3},\mathbf {a} _{3}\rangle &\ldots \\\vdots &\vdots &\vdots &\ddots \end{pmatrix}}.}$

Note que ${\displaystyle \langle \mathbf {e} _{j},\mathbf {a} _{j}\rangle =\|\mathbf {u} _{j}\|,}$ ${\displaystyle \langle \mathbf {e} _{j},\mathbf {a} _{k}\rangle =0\mathrm {~~para~~} j>k,}$ y ${\displaystyle QQ^{T}=I}$, entonces ${\displaystyle Q^{T}=Q^{-1}}$.

Si se considera la descomposición de

${\displaystyle A={\begin{pmatrix}12&-51&4\\6&167&-68\\-4&24&-41\end{pmatrix}}.}$

Se busca la matriz ortogonal ${\displaystyle Q}$ tal que

${\displaystyle {\begin{matrix}Q\,Q^{T}=I.\end{matrix}}}$

Por lo que calculamos ${\displaystyle Q}$ mediante Gram-Schmidt como sigue:

${\displaystyle U={\begin{pmatrix}\mathbf {u} _{1}&\mathbf {u} _{2}&\mathbf {u} _{3}\end{pmatrix}}={\begin{pmatrix}12&-69&-58/5\\6&158&6/5\\-4&30&-33\end{pmatrix}};}$

${\displaystyle Q={\begin{pmatrix}{\frac {\mathbf {u} _{1}}{\|\mathbf {u} _{1}\|}}&{\frac {\mathbf {u} _{2}}{\|\mathbf {u} _{2}\|}}&{\frac {\mathbf {u} _{3}}{\|\mathbf {u} _{3}\|}}\end{pmatrix}}={\begin{pmatrix}6/7&-69/175&-58/175\\3/7&158/175&6/175\\-2/7&6/35&-33/35\end{pmatrix}};}$

Por lo tanto, tenemos

${\displaystyle {\begin{matrix}A=Q\,Q^{T}A=QR;\end{matrix}}}$

${\displaystyle {\begin{matrix}R=Q^{T}A=\end{matrix}}{\begin{pmatrix}14&21&-14\\0&175&-70\\0&0&35\end{pmatrix}}.}$

Una transformación de Householder o reflexión de Householder es una transformación que refleja el espacio con respecto a un plano determinado. Esta propiedad se puede utilizar para realizar la transformación QR de una matriz si tenemos en cuenta que es posible elegir la matriz de Householder de manera que un vector elegido quede con una única componente no nula tras ser transformado (es decir, premultiplicando por la matriz de Householder). Gráficamente, esto significa que es posible reflejar el vector elegido respecto de un plano de forma que el reflejo quede sobre uno de los ejes de la base cartesiana.

La manera de elegir el plano de reflexión y formar la matriz de Householder asociada es el siguiente:

Sea ${\displaystyle \mathbf {x} }$ un vector columna arbitrario m-dimensional tal que ||${\displaystyle \mathbf {x} }$|| = |α|, donde α es un escalar; (si el algoritmo se implementa utilizando aritmética de coma flotante, entonces α debe adoptar el signo contrario que ${\displaystyle x}$₁ para evitar pérdida de precisión).

Entonces, siendo ${\displaystyle \mathbf {e} _{1}}$ el vector (1,0,...,0)^T, y ||·|| la norma euclídea, se define:

${\displaystyle \mathbf {u} =\mathbf {x} -\alpha \mathbf {e} _{1},}$

${\displaystyle \mathbf {v} ={\mathbf {u} \over ||\mathbf {u} ||},}$

${\displaystyle Q=I-2\mathbf {v} \mathbf {v} ^{T}.}$

${\displaystyle v}$ es un vector unitario perpendicular al plano de reflexión elegido. ${\displaystyle Q}$ es una matriz de Householder asociada a dicho plano.

${\displaystyle Qx=(\alpha ,0,\cdots ,0)^{T}.\,}$

Esta propiedad se puede usar para transformar gradualmente los vectores columna de una matriz A de dimensiones m por n en una matriz triangular superior. En primer lugar se multiplica A con la matriz de Householder Q₁ que obtenemos al elegir como vector ${\displaystyle \mathbf {x} }$ la primera columna de la matriz. Esto proporciona una matriz QA con ceros en la primera columna (excepto el elemento de la primera fila).

${\displaystyle Q_{1}A={\begin{bmatrix}\alpha _{1}&\star &\dots &\star \\0&&&\\\vdots &&A'&\\0&&&\end{bmatrix}}}$

el procedimiento se puede repetir para A′ (que se obtiene de A eliminando la primera fila y columna), obteniendo así una matriz de Householder Q′₂. Hay que tener en cuenta que Q′₂ es menor que Q₁. Para conseguir que esta matriz opere con Q₁A en lugar de A′ es necesario expandirla hacia arriba a la izquierda, completando con un uno en la diagonal, o en general:

${\displaystyle Q_{k}={\begin{pmatrix}I_{k-1}&0\\0&Q_{k}'\end{pmatrix}}.}$

Tras repetir el proceso ${\displaystyle t}$ veces, donde ${\displaystyle t=\min(m-1,n)}$,

${\displaystyle R=Q_{t}\cdots Q_{2}Q_{1}A}$

es una matriz triangular superior. De forma que tomando

${\displaystyle Q=Q_{1}Q_{2}\cdots Q_{t}}$

${\displaystyle A=Q^{T}R}$ es una descomposición QR de la matriz ${\displaystyle A}$.

Este método tiene una estabilidad numérica mayor que la del método de Gram-Schmidt descrito arriba.

Una pequeña variación de este método se utiliza para obtener matrices semejantes con la forma de Hessenberg, muy útiles en el cálculo de autovalores por acelerar la convergencia del algoritmo QR reduciendo así enormemente su coste computacional.

Vamos a calcular la descomposición de la matriz

${\displaystyle A={\begin{pmatrix}12&-51&4\\6&167&-68\\-4&24&-41\end{pmatrix}}.}$

En primer lugar necesitamos encontrar una reflexión que transforme la primera columna de la matriz A, vector ${\displaystyle \mathbf {a} _{1}=(12,6,-4)^{T}}$, en ${\displaystyle \|\mathbf {a} _{1}\|\;\mathrm {e} _{1}=(14,0,0)^{T}.}$

usando la expresión,

${\displaystyle \mathbf {u} =\mathbf {a} _{1}-\alpha \mathbf {e} _{1},}$

y

${\displaystyle \mathbf {v} ={\mathbf {u} \over |\mathbf {u} |},}$

en nuestro caso :

${\displaystyle \alpha =14}$ y ${\displaystyle \mathbf {a} _{1}=(12,6,-4)^{T}}$

Por lo tanto

${\displaystyle \mathbf {u} =(-2,6,-4)^{T}}$ y ${\displaystyle \mathbf {v} ={1 \over {\sqrt {14}}}(-1,3,-2)^{T}}$, entonces

${\displaystyle Q_{1}=I-{2 \over {\sqrt {14}}{\sqrt {14}}}{\begin{pmatrix}-1\\3\\-2\end{pmatrix}}{\begin{pmatrix}-1&3&-2\end{pmatrix}}}$

${\displaystyle =I-{1 \over 7}{\begin{pmatrix}1&-3&2\\-3&9&-6\\2&-6&4\end{pmatrix}}}$

${\displaystyle ={\begin{pmatrix}6/7&3/7&-2/7\\3/7&-2/7&6/7\\-2/7&6/7&3/7\\\end{pmatrix}}.}$

Ahora observamos:

${\displaystyle Q_{1}A={\begin{pmatrix}14&21&-14\\0&-49&-14\\0&168&-77\end{pmatrix}},}$

con lo que ya casi tenemos una matriz triangular. Sólo necesitamos hacer cero en el elemento (3,2).

Tomando la submatriz bajo el (1, 1) y aplicando de nuevo el proceso a

${\displaystyle A'=M_{11}={\begin{pmatrix}-49&-14\\168&-77\end{pmatrix}}.}$

Mediante el mismo método que antes obtenemos la matriz de Householder

${\displaystyle Q_{2}={\begin{pmatrix}1&0&0\\0&-7/25&24/25\\0&24/25&7/25\end{pmatrix}}}$

Finalmente obtenemos

${\displaystyle Q=Q_{1}Q_{2}={\begin{pmatrix}6/7&-69/175&58/175\\3/7&158/175&-6/175\\-2/7&6/35&33/35\end{pmatrix}}}$

${\displaystyle R=Q^{\top }A={\begin{pmatrix}14&21&-14\\0&175&-70\\0&0&35\end{pmatrix}}.}$

La matriz Q es ortogonal y R es triangular superior, de forma que A = QR es la descomposición QR buscada.

Las descomposiciones QR también puden calcularse utilizando una serie de rotaciones de Givens. Cada rotación anula (hace cero) un elemento en la subdiagonal de la matriz, formando de este modo la matriz R. La concatenación de todas las rotaciones de Givens realizadas, forma la matriz ortogonal Q.

En la práctica, las rotaciones de Givens no se utilizan en la actualidad para construir una matriz completa y realizar un producto de matrices. En su lugar, se utiliza un procedimiento de rotación de Givens, que es equivalente a la multiplicación reducida de matrices de Givens, sin el trabajo extra de manejar los elementos reducidos. El procedimiento de rotación de Givens es útil en situaciones donde sólo pocos elementos fuera de la diagonal necesitan ser anulados y es más fácil de paralelizar que las transformaciones de Householder.

Calculemos la descomposición de

${\displaystyle A={\begin{pmatrix}12&-51&4\\6&167&-68\\-4&24&-41\end{pmatrix}}}$

Primero, necesitamos formar una matriz de rotación tal que hagamos cero el elemento más inferior a la izquierda, ${\displaystyle \mathbf {a} _{31}=-4}$. Construimos esta matriz empleando el método de la rotación de Givens y llamamos la matriz resultante ${\displaystyle G_{1}}$. Rotamos primero el vector ${\displaystyle (6,-4)}$, representándolo a lo largo del eje X. Este vector forma un ángulo ${\displaystyle \theta =\arctan({-4 \over 6})}$. Creamos la matriz ortogonal de rotación de Givens, ${\displaystyle G_{1}}$:

${\displaystyle G_{1}={\begin{pmatrix}1&0&0\\0&\cos(\theta )&\sin(\theta )\\0&-\sin(\theta )&\cos(\theta )\end{pmatrix}}}$

${\displaystyle \approx {\begin{pmatrix}1&0&0\\0&0.83205&-0.55470\\0&0.55470&0.83205\end{pmatrix}}}$

Y el resultado de ${\displaystyle G_{1}A}$ tiene ahora un cero en el ${\displaystyle \mathbf {a} _{31}}$ elemento.

${\displaystyle G_{1}A\approx {\begin{pmatrix}12&-51&4\\7.21110&125.63959&-33.83671\\0&112.60414&-71.83368\end{pmatrix}}}$

Procedemos análogamente con las matrices de Givens ${\displaystyle G_{2}}$ y ${\displaystyle G_{3}}$, que hacen cero los elementos subdiagonales ${\displaystyle a_{21}}$ y ${\displaystyle a_{32}}$, formando una matriz triangular ${\displaystyle R}$. La matriz ortogonal ${\displaystyle Q^{T}}$ es formada a partir del producto en cadena de todas las matrices de Givens ${\displaystyle Q^{T}=G_{3}G_{2}G_{1}}$. Luego tenemos ${\displaystyle G_{3}G_{2}G_{1}A=Q^{T}A=R}$, y la descomposición QR es ${\displaystyle A=QR}$.