1. Introducción y contexto

El método de Gauss‑Seidel es una técnica iterativa clásica para resolver sistemas lineales de la forma Ax = b. Se emplea cuando la matriz A es grande, densa o dispersa y los métodos directos (LU, Cholesky) resultan costosos en tiempo o memoria. Es especialmente útil en problemas de simulación numérica, optimización y modelado de procesos físicos (por ejemplo, ecuaciones de difusión).

Este artículo cubre:

• Fundamentos matemáticos y condiciones de convergencia.
• Implementación paso a paso en Python 3 (sin dependencias externas y con NumPy).
• Comparativa con el método de Jacobi.
• Optimización para matrices dispersas y paralelismo.
• Diagnóstico de errores comunes y mejores prácticas.

2. Fundamento teórico

Partimos de la descomposición de A en sus componentes diagonal (D), estrictamente inferior (L) y estrictamente superior (U):

A = D + L + U

El método de Gauss‑Seidel actualiza la solución x usando la fórmula:

x^{(k+1)} = (D + L)^{-1}\bigl(b - U x^{(k)}\bigr)

En forma component‑wise, para cada fila i:

x_i^{(k+1)} = \frac{1}{a_{ii}}\Bigl(b_i - \sum_{j=1}^{i-1} a_{ij}x_j^{(k+1)} - \sum_{j=i+1}^{n} a_{ij}x_j^{(k)}\Bigr)

Observa que los valores ya actualizados del mismo ciclo (k+1) se usan inmediatamente, lo que suele acelerar la convergencia respecto a Jacobi.

Condiciones de convergencia (criterio suficiente):

• Matriz estrictamente diagonal dominante: |a_{ii}| > \sum_{j\neq i}|a_{ij}| para todo i.
• Matriz simétrica y definida positiva (SPD). En este caso, Gauss‑Seidel está garantizado a converger.

3. Implementación en Python

A continuación se muestra una versión mínima, clara y altamente legible del algoritmo.

import numpy as np
def gauss_seidel(A, b, x0=None, tol=1e-8, max_iter=500):
    """Resuelve Ax = b usando Gauss‑Seidel.
    Parámetros
    ----------
    A : ndarray (n, n)
        Matriz coeficiente (debe ser cuadrada).
    b : ndarray (n,)
        Vector del lado derecho.
    x0 : ndarray (n,), opcional
        Aproximación inicial. Si es None se usa un vector de ceros.
    tol : float, opcional
        Tolerancia de parada basada en la norma L2 del residuo.
    max_iter : int, opcional
        Número máximo de iteraciones.
    Returns
    -------
    x : ndarray (n,)
        Aproximación de la solución.
    info : dict
        Información de convergencia (iteraciones, residuo final, éxito).
    """
    n = A.shape[0]
    if x0 is None:
        x = np.zeros_like(b, dtype=float)
    else:
        x = x0.astype(float)
    for k in range(1, max_iter + 1):
        x_old = x.copy()
        for i in range(n):
            sigma = np.dot(A[i, :i], x[:i]) + np.dot(A[i, i + 1 :], x_old[i + 1 :])
            x[i] = (b[i] - sigma) / A[i, i]
        # Cálculo del residuo
        residual = np.linalg.norm(b - A @ x)
        if residual < tol:
            return x, {"iterations": k, "residual": residual, "converged": True}
    return x, {"iterations": max_iter, "residual": residual, "converged": False}

Esta función está diseñada para ser compatible con matrices densas y dispersas (usando scipy.sparse se pueden pasar vistas de fila con .toarray() o .dot()).

4. Ejemplo práctico paso a paso

Resolvamos el siguiente sistema 3x3:

A = [[ 4, -1,  0],
     [-1,  4, -1],
     [ 0, -1,  3]]
b = [15, 10, 10]

La matriz es estrictamente diagonal dominante, por lo que Gauss‑Seidel convergerá rápidamente.

import numpy as np
A = np.array([[4, -1, 0],
              [-1, 4, -1],
              [0, -1, 3]], dtype=float)
b = np.array([15, 10, 10], dtype=float)
x, info = gauss_seidel(A, b, tol=1e-10)
print("Solución:", x)
print("Info:", info)

Salida esperada:

Solución: [4.99999999 5.00000001 5.        ]
Info: {'iterations': 13, 'residual': 3.2e-11, 'converged': True}

Obsérvese cómo la solución es [5, 5, 5] con una precisión de 1e-10 en solo 13 iteraciones.

5. Comparativa con el método de Jacobi

En la práctica, para sistemas moderados y con recursos de CPU limitados, Gauss‑Seidel suele ser la opción preferida. Cuando se dispone de arquitectura altamente paralela (GPU, clusters), Jacobi o variantes como Red‑Black Gauss‑Seidel pueden ser más adecuadas.

6. Optimización para matrices dispersas

Los problemas reales (por ejemplo, discretizaciones de PDE) generan matrices con miles o millones de filas pero con muy pocos valores no nulos por fila. Utilizar scipy.sparse reduce drásticamente la huella de memoria y mejora la velocidad.

from scipy.sparse import csr_matrix
# Construcción de una matriz tridiagonal 10000x10000
n = 10000
main = 2 * np.ones(n)
off  = -1 * np.ones(n-1)
A_sparse = csr_matrix(np.diag(main) + np.diag(off, k=1) + np.diag(off, k=-1))
b = np.ones(n)
# Versión optimizada del algoritmo
def gauss_seidel_sparse(A, b, x0=None, tol=1e-8, max_iter=2000):
    n = A.shape[0]
    if x0 is None:
        x = np.zeros(n)
    else:
        x = x0.copy()
    D_inv = 1.0 / A.diagonal()
    for k in range(max_iter):
        x_old = x.copy()
        for i in range(n):
            row_start = A.indptr[i]
            row_end   = A.indptr[i+1]
            Ai = A.indices[row_start:row_end]
            Av = A.data[row_start:row_end]
            sigma = 0.0
            for idx, j in enumerate(Ai):
                if j != i:
                    sigma += Av[idx] * x[j] if j < i else Av[idx] * x_old[j]
            x[i] = (b[i] - sigma) * D_inv[i]
        if np.linalg.norm(b - A @ x) < tol:
            break
    return x
x = gauss_seidel_sparse(A_sparse, b)
print("Residuo final:", np.linalg.norm(b - A_sparse @ x))

Esta variante evita la creación de matrices densas y mantiene la complejidad O(nnz), donde nnz es el número de elementos no nulos.

7. Depuración y troubleshooting

• División por cero: Asegúrese de que a_{ii} ≠ 0 para todas las filas. Si la diagonal contiene ceros, reordene el sistema o aplique pivoting parcial.
• Falta de convergencia: Verifique la condición de diagonal dominante o la positividad definida. En caso contrario, considere pre‑condicionar el sistema (por ejemplo, multiplicando por D^{-1}).
• Oscilaciones: Si la norma del residuo se estabiliza sin bajar, pruebe relajación sucesiva (SOR) con un factor ω entre 1 y 2.
• Rendimiento bajo: Use versiones vectorizadas de NumPy o numba JIT para acelerar los bucles internos.

8. Extensiones avanzadas

8.1. Método SOR (Successive Over‑Relaxation)

Introduce un factor de relajación ω que acelera la convergencia:

x_i^{(k+1)} = (1-ω) x_i^{(k)} + \frac{ω}{a_{ii}}\Bigl(b_i - \sum_{ji} a_{ij}x_j^{(k)}\Bigr)

El valor óptimo de ω depende de la condición espectral de A; típicamente se elige entre 1.0 y 1.9.

8.2. Red‑Black Gauss‑Seidel (para GPU)

Divide la malla en dos conjuntos (rojo y negro) de forma que cada punto rojo depende solo de negros y viceversa. Esto permite paralelizar cada “color” en GPU o en arquitecturas multi‑core.

9. Buenas prácticas y checklist final

1. Compruebe la diagonal dominante o SPD antes de lanzar el algoritmo.
2. Use float64 para evitar pérdida de precisión en problemas mal condicionados.
3. Establezca una tolerancia realista (1e-8 a 1e-12) según la escala de b.
4. Implemente un límite máximo de iteraciones y registre el número de iteraciones realizadas.
5. Para matrices grandes, prefiera la variante scipy.sparse o una implementación con numba para el bucle interno.
6. Si la convergencia es lenta, pruebe SOR o cambie a un solver directo (e.g., scipy.sparse.linalg.spsolve) como referencia.