Algoritmo de Descomposición QR
Todo lo que necesitas saber para entender, implementar y aplicar la factorización QR en Python.
1. Fundamento Matemático
La descomposición QR factoriza una matriz A (de dimensión m × n, m ≥ n) en el producto:
A = Q·R
- Q: matriz ortogonal (QᵀQ = I) de tamaño m × n (o m × m si se completa).
- R: matriz triangular superior de tamaño n × n.
Esta factorización es la base de algoritmos de resolución de sistemas lineales, cálculo de valores propios y regresión lineal.
2. Algoritmos Clásicos
Existen varios métodos para obtener Q y R:
- Gram‑Schmidt clásico (intuitivo pero numéricamente inestable).
- Gram‑Schmidt modificado (mejora la estabilidad).
- Householder reflections (el método por defecto en la mayoría de librerías).
- Givens rotations (ideal para matrices dispersas).
3. Implementación en Python
A continuación se presentan tres enfoques: usando numpy.linalg.qr, scipy.linalg.qr (con opciones de Householder y Givens) y una implementación manual con Gram‑Schmidt modificado.
3.1. Usando NumPy (Householder)
import numpy as np
# Matriz de ejemplo
A = np.array([[12, -51, 4],
[ 6, 167, -68],
[-4, 24, -41]], dtype=float)
Q, R = np.linalg.qr(A)
print("Q =\n", Q)
print("R =\n", R)
NumPy emplea reflecciones de Householder, ofreciendo alta precisión y rendimiento O(m·n²).
3.2. Usando SciPy (Elección de algoritmo)
import numpy as np, scipy.linalg as la
A = np.random.rand(5, 3)
# QR con Householder (default)
Q_h, R_h = la.qr(A, mode='economic')
# QR con Givens (más eficiente para matrices dispersas)
Q_g, R_g = la.qr(A, mode='economic', pivoting=False, check_finite=False, overwrite_a=True)
print("Q (Householder) shape:", Q_h.shape)
print("Q (Givens) shape:", Q_g.shape)
SciPy permite especificar mode='economic' para obtener Q de tamaño m×n y R de n×n.
3.3. Implementación Manual – Gram‑Schmidt Modificado
import numpy as np
def gram_schmidt_mod(A):
"""Devuelve Q y R usando Gram‑Schmidt modificado.
A: matriz (m, n) con m >= n.
"""
m, n = A.shape
Q = np.zeros((m, n))
R = np.zeros((n, n))
for k in range(n):
v = A[:, k]
for j in range(k):
R[j, k] = np.dot(Q[:, j].T, A[:, k])
v = v - R[j, k] * Q[:, j]
R[k, k] = np.linalg.norm(v)
if R[k, k] == 0:
raise ValueError('Matriz singular o columnas linealmente dependientes')
Q[:, k] = v / R[k, k]
return Q, R
A = np.array([[1., 1., 0.], [1., 0., 1.], [0., 1., 1.]])
Q, R = gram_schmidt_mod(A)
print('Q =\n', Q)
print('R =\n', R)
print('Verificación A ≈ Q·R: \n', np.allclose(A, Q @ R))
Esta versión es didáctica pero menos estable que Householder; se recomienda solo para educación.
4. Comparativa con Otras Factorizaciones
| Método | Tipo de salida | Complejidad | Estabilidad Numérica | Casos de uso típicos |
|---|---|---|---|---|
| QR (Householder) | Q ortogonal, R triangular | O(m·n²) | Alta | Resolución de LS, eigen‑decomposition, ortogonalización |
| QR (Givens) | Q ortogonal, R triangular | O(k·n) para matrices dispersas (k = número de rotaciones) | Alta | Matrices muy dispersas, hardware SIMD |
| LU | L triangular, U triangular | O(n³) (para cuadradas) | Media (requiere pivoteo) | Sistemas lineales directos, determinantes |
| SVD | U, Σ, Vᵀ | O(m·n²) (similar a QR) pero con mayor constante | Muy alta | Compresión, análisis de rango, pseudoinversas |
5. Buenas Prácticas y Optimización
- Preferir implementaciones de librerías (NumPy/SciPy) para aprovechar BLAS/LAPACK optimizado y paralelismo OpenBLAS/MKL.
- En entornos de gran escala, usar
scipy.linalg.qrconoverwrite_a=Truepara reducir uso de memoria. - Para matrices extremadamente dispersas (> 90% ceros), emplear
scipy.sparse.linalg.spsolvejunto a Givens o factorizaciones QR sparsas (scipy.sparse.linalg.qrexperimental). - Validar la ortogonalidad de Q con
np.allclose(Q.T @ Q, np.eye(Q.shape[1]))antes de usarla en algoritmos críticos. - En problemas de regresión, combinar QR con
numpy.linalg.lstsqpara evitar la formación explícita de la pseudo‑inversa.
6. Solución de Problemas Comunes (Troubleshooting)
- Q no es ortogonal
- Posibles causas: uso de algoritmo Gram‑Schmidt sin re‑ortogonalización o datos mal condicionados. Solución: cambiar a
numpy.linalg.qr(Householder) o aplicar re‑ortogonalización conscipy.linalg.orth. - R tiene ceros en la diagonal
- Indica columnas linealmente dependientes. Eliminar columnas redundantes o usar
qrcon pivoteo (mode='full') para obtener permutación que mejore la diagonal. - Desbordamiento de memoria
- Trabajar con matrices muy grandes (> 10⁴×10⁴). Solución: usar versiones esparsas o dividir la operación en bloques (algoritmo de bloque Householder).
© 2025 BlogTech – Todos los derechos reservados.
Algoritmo de Descomposición QR: Teoría, Implementación y Ejemplos en Python