Introducción

Una base ortonormal es un conjunto de vectores que son unitarios y mutuamente ortogonales. Estas bases son esenciales en áreas como procesamiento de señales, machine learning, optimización y métodos numéricos porque simplifican cálculos de proyección, descomposición y reducción dimensional.

En este artículo veremos:

• Fundamentos teóricos del algoritmo de Gram‑Schmidt y sus variantes.
• Implementaciones en NumPy y SciPy.
• Comparativa con QR y SVD.
• Casos de uso reales y buenas prácticas de rendimiento y precisión.

Fundamentos teóricos

Dados n vectores linealmente independientes v₁, v₂, …, vₙ ∈ ℝᵐ, el proceso de Gram‑Schmidt genera una base ortonormal u₁, u₂, …, uₙ mediante:

u₁ = v₁ / ‖v₁‖
for k = 2 … n:
    wₖ = vₖ - Σ_{j=1}^{k-1} (⟨vₖ, uⱼ⟩) uⱼ
    uₖ = wₖ / ‖wₖ‖

Donde ⟨·,·⟩ es el producto interno y ‖·‖ la norma euclídea.

La versión modificada (MGS) mejora la estabilidad numérica al actualizar vₖ en cada iteración, reduciendo la pérdida de ortogonalidad en precisión finita.

Implementación en Python

import numpy as np
def gram_schmidt(A):
    """Devuelve una base ortonormal a partir de la matriz A (m x n)."""
    m, n = A.shape
    Q = np.zeros((m, n))
    for k in range(n):
        v = A[:, k]
        for j in range(k):
            proj = np.dot(Q[:, j], v) * Q[:, j]
            v = v - proj
        norm = np.linalg.norm(v)
        if norm < 1e-12:
            raise ValueError('Los vectores son linealmente dependientes')
        Q[:, k] = v / norm
    return Q
# Ejemplo de uso
A = np.random.rand(5, 3)
Q = gram_schmidt(A)
print('Q.T @ Q =\n', np.round(Q.T @ Q, 12))  # debe ser ~I

def modified_gram_schmidt(A):
    m, n = A.shape
    Q = A.astype(float).copy()
    for i in range(n):
        norm = np.linalg.norm(Q[:, i])
        if norm < 1e-12:
            raise ValueError('Dependencia lineal detectada')
        Q[:, i] = Q[:, i] / norm
        for j in range(i+1, n):
            proj = np.dot(Q[:, i], Q[:, j])
            Q[:, j] = Q[:, j] - proj * Q[:, i]
    return Q
Qm = modified_gram_schmidt(A)
print('Ortogonalidad MGS:', np.max(np.abs(Qm.T @ Qm - np.eye(n))))

3️⃣ Usando scipy.linalg.qr (QR de Householder)

from scipy.linalg import qr
Q_qr, R = qr(A, mode='economic')
print('Ortogonalidad QR:', np.max(np.abs(Q_qr.T @ Q_qr - np.eye(n))))

Este método es el estándar de facto: rápido, estable y con soporte para float32, float64 e incluso tipos complejos.

Comparativa de rendimiento y precisión

Casos de uso en el mundo real

Buenas prácticas y consideraciones de seguridad numérica

• Escalado previo: Normaliza los datos a rango [0,1] o [-1,1] antes de aplicar Gram‑Schmidt para evitar overflow/underflow.
• Tipo de dato: Usa float64 en cálculos críticos; float32 puede ser suficiente en inferencia de IA cuando se combina con torch.nn.utils.prune.
• Re‑ortogonalización periódica: En algoritmos iterativos (p.ej., método de potencias) reaplica QR cada k iteraciones para controlar la pérdida de ortogonalidad.
• Detección de degeneración: Verifica que la norma de los vectores intermedios no caiga bajo un umbral (1e‑12) y, en caso contrario, aplica pivotado o elimina columnas dependientes.
• Seguridad de la memoria: Cuando se trabaja con datos sensibles (ej. biométricos), sobrescribe buffers temporales con numpy.fill(0) después de su uso para evitar filtraciones.

Optimización avanzada

Para matrices de gran escala (m,n > 10⁴) se recomiendan los siguientes enfoques:

1. BLAS/LAPACK multihilo: Usa numpy.linalg.qr con OpenBLAS o Intel MKL compilados para aprovechar todos los núcleos.
2. Algoritmos aleatorios: Randomized QR (disponible en scikit‑learn.utils.extmath.randomized_range_finder) reduce la complejidad a O(m·n·log(k)) para obtener una base de rango k.
3. GPU: Bibliotecas como cupy.linalg.qr o torch.linalg.qr trasladan la carga a la GPU, logrando mejoras de 5‑10× en tiempo.
4. Batch processing: Cuando se procesan cientos de matrices pequeñas, agrúpalas en un tensor 3‑D y usa torch.linalg.qr con batch_first=True para paralelismo.

Troubleshooting frecuente