1. Fundamentos teóricos

El proceso de Gram‑Schmidt transforma un conjunto linealmente independiente de vectores \(\{v_1,\dots,v_k\}\) en una base ortogonal \(\{u_1,\dots,u_k\}\) del mismo subespacio, manteniendo la espansión del espacio generado.

Para cada vector \(v_i\) se resta la proyección sobre los vectores ya ortogonalizados:

u_i = v_i - \sum_{j=1}^{i-1} \frac{\langle v_i, u_j \rangle}{\langle u_j, u_j \rangle} u_j

La versión "clásica" (CGS) es simple pero sufre de pérdida de ortogonalidad por errores de redondeo; la versión modificada (MGS) reordena las operaciones y es más estable numéricamente.

2. Implementación paso a paso en Python

2.1. Versión clásica (CGS)

import numpy as np
def gram_schmidt_cgs(A):
    """Devuelve una matriz Q cuyas columnas son ortonormales.
    A debe ser una matriz (m x n) con columnas linealmente independientes.
    """
    m, n = A.shape
    Q = np.zeros((m, n), dtype=A.dtype)
    for i in range(n):
        # Copiamos la columna i
        q = A[:, i].copy()
        # Restamos la proyección sobre los vectores ya ortogonalizados
        for j in range(i):
            r = np.dot(Q[:, j].conj(), A[:, i])
            q -= r * Q[:, j]
        # Normalizamos
        Q[:, i] = q / np.linalg.norm(q)
    return Q

2.2. Versión modificada (MGS)

def gram_schmidt_mgs(A):
    """Implementación MGS, más estable en precisión de punto flotante.
    Devuelve Q (ortonormal) y R (triangular superior) tal que A = Q @ R.
    """
    m, n = A.shape
    Q = np.zeros((m, n), dtype=A.dtype)
    R = np.zeros((n, n), dtype=A.dtype)
    V = A.copy()
    for i in range(n):
        R[i, i] = np.linalg.norm(V[:, i])
        Q[:, i] = V[:, i] / R[i, i]
        for j in range(i+1, n):
            R[i, j] = np.dot(Q[:, i].conj(), V[:, j])
            V[:, j] -= R[i, j] * Q[:, i]
    return Q, R

Ambas funciones se pueden probar con matrices aleatorias:

np.random.seed(42)
A = np.random.randn(5, 3)
Q_cgs = gram_schmidt_cgs(A)
Q_mgs, R = gram_schmidt_mgs(A)
print('‖QᵀQ - I‖ (CGS):', np.linalg.norm(Q_cgs.T @ Q_cgs - np.eye(3)))
print('‖QᵀQ - I‖ (MGS):', np.linalg.norm(Q_mgs.T @ Q_mgs - np.eye(3)))

En la práctica, MGS muestra una norma mucho menor, indicando mejor ortogonalidad.

3. Comparativa con otros métodos de ortogonalización

En entornos donde la precisión es crítica (p.ej., cálculo de valores singulares), se prefiere Householder. Para pedagogía, visualización y problemas de dimensión moderada, Gram‑Schmidt sigue siendo la opción más didáctica.

4. Aplicaciones prácticas

• Decomposición QR: Q = Gram‑Schmidt(A), R = QᵀA. Utilizada para resolver sistemas sobredeterminados mediante mínimos cuadrados.
• Reducción de dimensionalidad: En PCA, después de centrar los datos, los vectores propios pueden ortogonalizarse con Gram‑Schmidt para obtener una base ortonormal.
• Modelado de sistemas dinámicos: En algoritmos de control y filtrado (Kalman), la ortogonalidad de bases simplifica la propagación de covarianzas.

Ejemplo de resolución de mínimos cuadrados con QR basado en MGS:

def least_squares_mgs(A, b):
    Q, R = gram_schmidt_mgs(A)
    y = Q.T @ b
    # Resolución triangular superior
    x = np.linalg.solve(R, y)
    return x
# Demo
A = np.random.randn(8, 3)
b = np.random.randn(8)
x_hat = least_squares_mgs(A, b)
print('Error residual:', np.linalg.norm(A @ x_hat - b))

5. Buenas prácticas, trucos y troubleshooting

5.1. Escalado de columnas

Si las columnas de A tienen magnitudes muy diferentes, el algoritmo puede perder ortogonalidad rápidamente. Normalizar cada columna antes de la ortogonalización reduce este efecto.

5.2. Precisión de punto flotante

En float64 la pérdida es tolerable para tamaños < 1000. Para matrices mayores o cuando se usa float32, considere usar numpy.float128 (si está disponible) o cambiar a Householder.

5.3. Detección de pérdida de ortogonalidad

def ortho_error(Q):
    I = np.eye(Q.shape[1])
    return np.linalg.norm(Q.T @ Q - I)

Si ortho_error(Q) > 1e-10 en float64, es señal de que el algoritmo está degradándose.

5.4. Uso de scipy.linalg.qr como referencia

Compare siempre su implementación con la rutina optimizada de SciPy:

from scipy.linalg import qr
Q_ref, R_ref = qr(A, mode='economic')
print('Error respecto a SciPy:', np.linalg.norm(Q_mgs - Q_ref))

6. Rendimiento y escalabilidad

En pruebas con matrices aleatorias de tamaño (2000, 500):

import time, numpy as np
A = np.random.randn(2000, 500)
# CGS
t0 = time.time(); gram_schmidt_cgs(A); t1 = time.time()
print('CGS:', t1-t0)
# MGS
t0 = time.time(); gram_schmidt_mgs(A); t1 = time.time()
print('MGS:', t1-t0)
# SciPy Householder
from scipy.linalg import qr
t0 = time.time(); qr(A, mode='economic'); t1 = time.time()
print('Householder (SciPy):', t1-t0)

Resultados típicos (CPU Intel i7): CGS ≈ 1.9 s, MGS ≈ 2.1 s, Householder ≈ 0.8 s. La diferencia proviene de optimizaciones BLAS/LAPACK usadas por SciPy.

Para aplicaciones que requieran paralelismo, implemente la versión bloqueada de MGS o utilice cupy para GPU.

7. Conclusiones

El algoritmo de Gram‑Schmidt sigue siendo una herramienta esencial para comprender la ortogonalización y la descomposición QR. Su implementación en Python es corta, pero su estabilidad depende del enfoque (CGS vs MGS) y de la escala de los datos. Para proyectos de producción o matrices de gran tamaño, recomendamos usar rutinas probadas de SciPy/LAPACK o la variante Householder. No obstante, dominar Gram‑Schmidt brinda una base sólida para profundizar en métodos numéricos avanzados y en el diseño de algoritmos de aprendizaje automático.