Algoritmos para Matrices Banda y Tridiagonales

En este artículo profundizamos en la teoría y práctica de dos tipos de matrices dispersas muy usadas en simulaciones numéricas: matrices banda y matrices tridiagonales. Incluimos algoritmos clásicos (LU para matrices banda y el método de Thomas para tridiagonales), ejemplos completos en Python y una tabla comparativa con sus principales características.

1. Conceptos Básicos

1.1 Matriz Banda

Una matriz banda tiene todos los elementos no nulos concentrados en un número limitado de diagonales alrededor de la diagonal principal. El ancho de banda se define por dos valores:

• l: número de sub‑diagonales (por debajo de la diagonal principal).
• u: número de super‑diagonales (por encima de la diagonal principal).

Ejemplo (l = 1, u = 2):

[[a00, a01, a02, 0  , 0  ],
 [a10, a11, a12, a13, 0  ],
 [0  , a21, a22, a23, a24],
 [0  , 0  , a32, a33, a34],
 [0  , 0  , 0  , a43, a44]]

1.2 Matriz Tridiagonal

Una matriz tridiagonal es un caso particular de matriz banda con l = u = 1. Sólo la diagonal principal y las dos adyacentes pueden contener valores distintos de cero.

[[b0, c0, 0 , 0 , …],
 [a1, b1, c1, 0 , …],
 [0 , a2, b2, c2, …],
 …]

2. Algoritmos Clásicos

2.1 Descomposición LU para Matrices Banda

La descomposición LU aprovecha la estructura de banda para reducir la complejidad de O(n³) a O(n·(l+u)²). El algoritmo procede fila a fila, manteniendo sólo los elementos dentro del ancho de banda.

def lu_band(A, l, u):
    n = A.shape[0]
    L = np.eye(n)
    U = np.zeros_like(A)
    for i in range(n):
        # U diagonal y super‑diagonales
        for j in range(i, min(i+u+1, n)):
            sum_ = sum(L[i, k] * U[k, j] for k in range(max(0, i-l), i))
            U[i, j] = A[i, j] - sum_
        # L sub‑diagonales
        for j in range(i+1, min(i+l+1, n)):
            sum_ = sum(L[j, k] * U[k, i] for k in range(max(0, i-l), i))
            L[j, i] = (A[j, i] - sum_) / U[i, i]
    return L, U

2.2 Método de Thomas (Tri‑Diagonal Matrix Algorithm - TDMA)

El método de Thomas es una versión optimizada del algoritmo de eliminación Gauss‑Jordan para matrices tridiagonales, con complejidad lineal O(n). Se divide en dos fases: forward elimination y back substitution.

def thomas(a, b, c, d):
    """Resuelve Ax = d donde A es tridiagonal.
    a: sub‑diagonal (n‑1 elementos)
    b: diagonal principal (n elementos)
    c: super‑diagonal (n‑1 elementos)
    d: vector del lado derecho (n elementos)
    """
    n = len(b)
    # Forward sweep
    c_ = np.empty(n-1)
    d_ = np.empty(n)
    c_[0] = c[0] / b[0]
    d_[0] = d[0] / b[0]
    for i in range(1, n-1):
        denom = b[i] - a[i-1] * c_[i-1]
        c_[i] = c[i] / denom
        d_[i] = (d[i] - a[i-1] * d_[i-1]) / denom
    d_[n-1] = (d[n-1] - a[n-2] * d_[n-2]) / (b[n-1] - a[n-2] * c_[n-2])
    # Back substitution
    x = np.empty(n)
    x[-1] = d_[n-1]
    for i in range(n-2, -1, -1):
        x[i] = d_[i] - c_[i] * x[i+1]
    return x

3. Comparativa de Características

4. Buenas Prácticas y Optimización

• Usar tipos de dato float64 o float32 según la precisión requerida.
• Evitar bucles Python puros en favor de numpy o numba para acelerar cálculos.
• Para matrices muy grandes (>10⁶), considerar scipy.sparse con formato dia_matrix o bsr_matrix.
• Validar la condición de diagonal dominante antes de aplicar Thomas; en caso contrario, usar pivoting o LU general.
• Perfilado con cProfile o line_profiler para identificar cuellos de botella.

5. Troubleshooting Común

6. Implementación Completa – Caso de Estudio

Resolución de la ecuación de calor 1‑D usando diferencias finitas (tridiagonal) y comparación con una matriz banda de ancho 3.

import numpy as np
from scipy.sparse import diags
# Parámetros
n = 1000          # número de nodos internos
L = 1.0           # longitud del dominio
alpha = 0.01      # coeficiente de difusión
dx = L/(n+1)
# Matriz tridiagonal (Thomas)
a = np.full(n-1, -alpha/dx**2)   # sub‑diagonal
b = np.full(n,  2*alpha/dx**2)    # diagonal principal
c = np.full(n-1, -alpha/dx**2)   # super‑diagonal
# Condiciones de frontera Dirichlet: u(0)=u(L)=0
f = np.zeros(n)                  # vector fuente (cero en este caso)
# Solución usando Thomas
u_tridiag = thomas(a, b, c, f)
# Matriz banda usando SciPy (ancho = 2 sub‑ y 2 super‑diagonales)
offsets = [-2, -1, 0, 1, 2]
values = [np.full(n, -alpha/dx**2) if off != 0 else np.full(n, 2*alpha/dx**2) for off in offsets]
A_band = diags(values, offsets, shape=(n, n)).toarray()
# Solución usando LU banda (función lu_band definida antes)
L_mat, U_mat = lu_band(A_band, l=2, u=2)
# Forward substitution L y = f
y = np.linalg.solve(L_mat, f)
# Back substitution U x = y
u_band = np.linalg.solve(U_mat, y)
print('Error L2 entre soluciones:', np.linalg.norm(u_tridiag-u_band))

En la práctica, la solución tridiagonal es ≈ 5‑10 veces más rápida que la descomposición LU banda para n = 10⁶, mientras que el consumo de memoria es mucho menor al usar solo tres vectores.

7. Conclusiones

Las matrices banda y tridiagonales son pilares en la computación científica. Elegir el algoritmo adecuado (LU banda vs. Thomas) depende del ancho de banda, del tamaño del problema y de los requisitos de rendimiento. Con las técnicas y buenas prácticas presentadas, puedes construir soluciones robustas, escalables y altamente optimizadas.