Introducción

La eliminación gaussiana (también conocida como método de eliminación de Gauss) es la técnica fundamental para resolver sistemas de ecuaciones lineales de la forma Ax = b. Consiste en transformar la matriz A en una forma triangular superior mediante operaciones elementales de fila, y luego aplicar sustitución regresiva.

Este artículo muestra cómo implementarla en Python desde cero, con NumPy y SciPy, y compara su rendimiento y estabilidad con métodos alternativos como LU y QR.

Fundamentos Teóricos

• Pivoteo parcial: intercambio de filas para evitar divisiones por números cercanos a cero.
• Pivoteo total (opcional): busca el mayor valor absoluto en la submatriz restante.
• Complejidad típica: O(n³) para una matriz n×n.

La versión más segura en práctica incluye pivoteo parcial, pues reduce el error de redondeo y evita division by zero.

1. Implementación en Python puro

Esta versión sirve para entender cada paso del algoritmo sin depender de librerías externas.

def gauss_elimination(A, b):
    """Resuelve Ax = b usando eliminación gaussiana con pivoteo parcial.
    Parámetros:
        A (list[list[float]]): matriz de coeficientes (copia será modificada).
        b (list[float]): vector de términos independientes.
    Retorna:
        x (list[float]): solución del sistema.
    """
    n = len(A)
    # --- Paso 1: Transformar a triangular superior ---
    for k in range(n):
        # Pivoteo parcial
        max_row = max(range(k, n), key=lambda i: abs(A[i][k]))
        if A[max_row][k] == 0:
            raise ValueError("Matriz singular o sistema sin solución única")
        # Intercambio de filas
        A[k], A[max_row] = A[max_row], A[k]
        b[k], b[max_row] = b[max_row], b[k]
        # Eliminación de filas inferiores
        for i in range(k + 1, n):
            factor = A[i][k] / A[k][k]
            for j in range(k, n):
                A[i][j] -= factor * A[k][j]
            b[i] -= factor * b[k]
    # --- Paso 2: Sustitución regresiva ---
    x = [0.0] * n
    for i in reversed(range(n)):
        sum_ax = sum(A[i][j] * x[j] for j in range(i + 1, n))
        x[i] = (b[i] - sum_ax) / A[i][i]
    return x
# Ejemplo de uso
A = [[2, 1, -1],
     [-3, -1, 2],
     [-2, 1, 2]]
b = [8, -11, -3]
solution = gauss_elimination([row[:] for row in A], b[:])
print(solution)  # → [2.0, 3.0, -1.0]

2. Implementación con NumPy

NumPy brinda operaciones vectorizadas que reducen el número de bucles y mejoran el rendimiento.

import numpy as np
def gauss_numpy(A, b):
    A = A.astype(float)
    b = b.astype(float)
    n = A.shape[0]
    for k in range(n):
        # Pivoteo parcial
        max_row = np.argmax(np.abs(A[k:, k])) + k
        if A[max_row, k] == 0:
            raise np.linalg.LinAlgError("Matriz singular")
        # Intercambio de filas
        if max_row != k:
            A[[k, max_row]] = A[[max_row, k]]
            b[[k, max_row]] = b[[max_row, k]]
        # Eliminación
        for i in range(k + 1, n):
            factor = A[i, k] / A[k, k]
            A[i, k:] -= factor * A[k, k:]
            b[i] -= factor * b[k]
    # Sustitución regresiva
    x = np.empty(n)
    for i in range(n - 1, -1, -1):
        x[i] = (b[i] - np.dot(A[i, i + 1:], x[i + 1:])) / A[i, i]
    return x
# Uso
A = np.array([[2, 1, -1],
              [-3, -1, 2],
              [-2, 1, 2]])
b = np.array([8, -11, -3])
print(gauss_numpy(A, b))  # → [ 2.  3. -1.]

Ventaja: Operaciones en bloque (A[i, k:]) aprovechan SIMD bajo el capó.

3. Uso de SciPy – Factorización LU (alternativa)

Si tu objetivo es rendimiento y robustez, la factorización LU de scipy.linalg es la opción recomendada. Internamente usa algoritmos LAPACK altamente optimizados.

from scipy.linalg import lu_factor, lu_solve
import numpy as np
A = np.array([[2, 1, -1],
              [-3, -1, 2],
              [-2, 1, 2]], dtype=float)
b = np.array([8, -11, -3], dtype=float)
lu, piv = lu_factor(A)          # Factoriza A = P·L·U
x = lu_solve((lu, piv), b)      # Resuelve el sistema
print(x)  # → [ 2.  3. -1.]

Esta solución es más rápida que la implementación manual para matrices de tamaño medio‑grande (n > 500) y ofrece mejor manejo de singularidad.

Comparativa de Algoritmos

Rendimiento y Benchmarking

El siguiente script compara los tiempos de ejecución para n = 800 usando las tres técnicas.

import numpy as np, time, scipy.linalg as la
n = 800
A = np.random.rand(n, n)
b = np.random.rand(n)
# 1. Gauss puro (NumPy) – solo para referencia (muy lento)
start = time.time()
gauss_numpy(A.copy(), b.copy())
print('Gauss (NumPy):', time.time() - start)
# 2. LU de SciPy
start = time.time()
lu, piv = la.lu_factor(A)
lu_solve((lu, piv), b)
print('LU (SciPy):', time.time() - start)
# 3. QR (Householder) – usando la función de SciPy
start = time.time()
Q, R = la.qr(A, mode='economic')
x = la.solve_triangular(R, Q.T @ b)
print('QR (SciPy):', time.time() - start)

En máquinas modernas (Intel i7‑12700K, 32 GB RAM) los resultados típicos son:

Como se observa, la factorización LU es la más rápida para sistemas cuadrados densos.

Buenas Prácticas y Solución de Problemas

• Escalado de la matriz: antes de aplicar pivoteo, escalar filas/columnas reduce la pérdida de precisión.
• Control de singularidad: usar np.linalg.cond(A) para detectar matrices mal condicionadas (cond > 1e12 sugiere problemas).
• Precisión de punto flotante: en problemas críticos, emplear numpy.float64 o mpmath para precisión arbitraria.
• Iteraciones de verificación: después de obtener x, comprobar np.allclose(A @ x, b) para validar la solución.
• Uso de bibliotecas optimizadas: BLAS/LAPACK (OpenBLAS, Intel MKL) aceleran operaciones de matrices en NumPy y SciPy.

Casos de Uso en el Mundo Real

La eliminación gaussiana y sus variantes aparecen en numerosos dominios:

• Ingeniería estructural: cálculo de fuerzas en sistemas de barra‑nudo.
• Simulación de circuitos eléctricos: análisis nodal (MNA) utiliza LU para resolver G·V = I.
• Economía y finanzas: modelos de regresión lineal múltiple (mínimos cuadrados) se resuelven con QR.
• Machine Learning: entrenamiento de regresión lineal y ridge regression.

Conclusión

La eliminación gaussiana sigue siendo una herramienta pedagógica esencial y, con pivoteo parcial, una solución viable para sistemas de tamaño moderado. Sin embargo, en entornos de producción y de alto rendimiento, la factorización LU (SciPy) o los métodos iterativos para matrices dispersas son generalmente preferibles.

Dominar estas técnicas y saber cuándo aplicar cada una te permitirá diseñar soluciones robustas, eficientes y escalables.