Sistemas de ecuaciones lineales con NumPy

Sistemas exactos con np.linalg.solve

Un sistema de ecuaciones lineales se expresa en forma matricial como $Ax = b$, donde $A$ es la matriz de coeficientes, $x$ el vector de incógnitas y $b$ el vector de términos independientes. Cuando $A$ es una matriz cuadrada con determinante distinto de cero, el sistema tiene una solución única.

La función np.linalg.solve resuelve este tipo de sistemas de forma eficiente utilizando factorización LU internamente, lo qué es numéricamente superior a calcular la inversa de la matriz.

import numpy as np

# Sistema de 2 ecuaciones:
#   3x + y = 9
#   x + 2y = 8
A = np.array([[3, 1],
              [1, 2]])
b = np.array([9, 8])

x = np.linalg.solve(A, b)
print(x)  # [2. 3.]
# Verificación: x=2, y=3 => 3*2+3=9, 2+2*3=8

Para verificar que la solución es correcta, se puede comprobar que $Ax$ es igual a $b$:

# Verificación numérica
residuo = A @ x - b
print(residuo)  # [0. 0.] (o valores muy cercanos a cero)
print(np.allclose(A @ x, b))  # True

Sistemas de mayor dimensión

El mismo procedimiento escala a sistemas con cualquier número de ecuaciones e incógnitas, siempre que la matriz sea cuadrada e invertible.

# Sistema de 3 ecuaciones con 3 incógnitas
#   2x + y - z  = 8
#   -3x - y + 2z = -11
#   -2x + y + 2z = -3
A = np.array([[ 2,  1, -1],
              [-3, -1,  2],
              [-2,  1,  2]])
b = np.array([8, -11, -3])

x = np.linalg.solve(A, b)
print(x)  # [2. 3. -1.]

La función np.linalg.solve es preferible a calcular la inversa ($x = A^{-1}b$) porque es más rápida y numéricamente estable. Calcular la inversa con np.linalg.inv duplica el coste computacional y amplifica los errores de redondeo.

Sistemas sobredeterminados con np.linalg.lstsq

En muchas aplicaciones prácticas el sistema tiene más ecuaciones que incógnitas (sistema sobredeterminado). En este caso no existe una solución exacta, pero se puede encontrar la solución que minimiza el error cuadrático, es decir, la solución por mínimos cuadrados.

La función np.linalg.lstsq resuelve el problema de minimización $\min_x |Ax - b|^2$.

import numpy as np

# 5 ecuaciones, 2 incógnitas (sobredeterminado)
A = np.array([[1, 1],
              [1, 2],
              [1, 3],
              [1, 4],
              [1, 5]])
b = np.array([2.1, 3.9, 6.2, 7.8, 10.1])

# rcond=None evita un warning de deprecación
x, residuos, rango, sv = np.linalg.lstsq(A, b, rcond=None)

print(f"Coeficientes: {x}")
print(f"Residuos: {residuos}")
print(f"Rango: {rango}")
print(f"Valores singulares: {sv}")

La función devuelve cuatro valores:

• x: la solución por mínimos cuadrados.
• residuos: la suma de los cuadrados de los residuos ($|Ax - b|^2$). Solo se devuelve cuando el sistema es sobredeterminado y la matriz tiene rango completo.
• rango: el rango efectivo de la matriz $A$.
• sv: los valores singulares de $A$.

flowchart TD
  subgraph Entrada[Tipo de sistema]
    E1["Sistema exacto<br>A cuadrada, det != 0"]
    E2["Sistema sobredeterminado<br>mas ecuaciones<br>que incognitas"]
    E3["Sistema subdeterminado<br>o singular"]
  end
  subgraph Función[Función recomendada]
    F1["np.linalg.solve"]
    F2["np.linalg.lstsq"]
    F3["np.linalg.pinv"]
  end
  E1 --> F1
  E2 --> F2
  E3 --> F3

Pseudoinversa con np.linalg.pinv

La pseudoinversa de Moore-Penrose ($A^+$) generaliza el concepto de inversa a matrices que no son cuadradas o que son singulares. Se calcula con np.linalg.pinv y permite resolver sistemas en los que np.linalg.solve no es aplicable.

import numpy as np

# Matriz no cuadrada (3x2)
A = np.array([[1, 2],
              [3, 4],
              [5, 6]])
b = np.array([1, 2, 3])

# Pseudoinversa
A_pinv = np.linalg.pinv(A)
print(f"Forma de A+: {A_pinv.shape}")  # (2, 3)

# Solución por pseudoinversa
x = A_pinv @ b
print(f"Solución: {x}")

La solución obtenida mediante la pseudoinversa coincide con la de np.linalg.lstsq para sistemas sobredeterminados. La pseudoinversa también es útil para matrices cuadradas singulares, donde np.linalg.solve lanzaría un error LinAlgError.

# Matriz singular (determinante = 0)
A_singular = np.array([[1, 2],
                        [2, 4]])

# np.linalg.solve fallaría aquí
# x = np.linalg.solve(A_singular, [3, 6])  # LinAlgError

# La pseudoinversa sí funciona
A_pinv = np.linalg.pinv(A_singular)
x = A_pinv @ np.array([3, 6])
print(f"Solución (norma mínima): {x}")

La pseudoinversa proporciona la solución de norma mínima cuando el sistema tiene infinitas soluciones, y la solución de mínimos cuadrados cuando no tiene solución exacta.

Número de condición y estabilidad numérica

El número de condición de una matriz indica cuánto se amplifica el error en la solución cuando hay pequeñas perturbaciones en los datos de entrada. Se calcula con np.linalg.cond.

import numpy as np

# Matriz bien condicionada
A_buena = np.array([[1, 0],
                     [0, 1]])
print(f"Cond(identidad): {np.linalg.cond(A_buena):.1f}")  # 1.0

# Matriz mal condicionada
A_mala = np.array([[1, 1],
                    [1, 1.0001]])
print(f"Cond(mal cond.): {np.linalg.cond(A_mala):.1f}")  # ~40000

Un número de condición cercano a 1 indica un sistema estable. Valores grandes (del orden de $10^{10}$ o superiores) indican que la solución puede ser poco fiable debido a la aritmética de punto flotante.

Efecto del mal condicionamiento

A = np.array([[1.0, 1.0],
              [1.0, 1.0001]])
b1 = np.array([2.0, 2.0001])
b2 = np.array([2.0, 2.0002])  # Perturbación mínima en b

x1 = np.linalg.solve(A, b1)
x2 = np.linalg.solve(A, b2)

print(f"x1: {x1}")  # [1. 1.]
print(f"x2: {x2}")  # [0. 2.] aprox. (cambio grande en x)
print(f"Diferencia en b: {np.linalg.norm(b2 - b1):.4f}")
print(f"Diferencia en x: {np.linalg.norm(x2 - x1):.4f}")

Una perturbación mínima en $b$ produce un cambio significativo en la solución, lo que confirma que el sistema está mal condicionado. Antes de confiar en la solución de un sistema lineal, conviene comprobar el número de condición.

cond = np.linalg.cond(A)
print(f"Número de condición: {cond:.1f}")

if cond > 1e10:
    print("Sistema mal condicionado: solución poco fiable")
else:
    print("Sistema bien condicionado")

Aplicación: regresión lineal por mínimos cuadrados

La regresión lineal se formula como un sistema sobredeterminado donde se buscan los coeficientes que mejor ajustan una recta (o un hiperplano) a un conjunto de puntos observados.

Dado un conjunto de puntos $(x_i, y_i)$, el modelo lineal $y = mx + c$ se puede expresar en forma matricial como:

$$ \begin{bmatrix} 1 & x_1 \ 1 & x_2 \ \vdots & \vdots \ 1 & x_n \end{bmatrix} \begin{bmatrix} c \ m \end{bmatrix} = \begin{bmatrix} y_1 \ y_2 \ \vdots \ y_n \end{bmatrix} $$

import numpy as np

# Datos: superficie (m2) y precio (miles de euros) de 8 pisos
superficie = np.array([50, 65, 80, 95, 110, 125, 140, 160])
precio = np.array([120, 155, 185, 220, 250, 280, 315, 360])

# Construir la matriz de diseño con columna de unos
A = np.column_stack([np.ones_like(superficie), superficie])

# Resolver por mínimos cuadrados
coefs, residuos, _, _ = np.linalg.lstsq(A, precio, rcond=None)
c, m = coefs

print(f"Ordenada: {c:.2f} miles EUR")
print(f"Pendiente: {m:.2f} miles EUR/m2")
print(f"Predicción para 100 m2: {c + m * 100:.2f} miles EUR")

Ajuste polinómico

Para relaciones no lineales se puede ajustar un polinomio de grado $n$ utilizando el mismo enfoque. La matriz de diseño se construye con las potencias de $x$.

# Datos con tendencia cuadrática
x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8])
y = np.array([1.0, 2.8, 7.1, 14.2, 23.5, 36.0, 50.8, 68.1, 88.5])

# Matriz de diseño para polinomio de grado 2: [1, x, x^2]
A = np.column_stack([np.ones_like(x), x, x**2])

coefs, _, _, _ = np.linalg.lstsq(A, y, rcond=None)
c0, c1, c2 = coefs

print(f"y = {c2:.3f}x^2 + {c1:.3f}x + {c0:.3f}")

# Predicciones
x_pred = np.array([9, 10])
y_pred = c0 + c1 * x_pred + c2 * x_pred**2
print(f"Predicciones: {np.round(y_pred, 1)}")

NumPy también ofrece np.polyfit y np.polynomial.polynomial.Polynomial.fit como atajos para el ajuste polinómico. Sin embargo, construir la matriz de diseño manualmente ofrece mayor flexibilidad y permite incluir términos personalizados como interacciones o funciones no polinómicas.

# Comparación con np.polyfit (devuelve coeficientes en orden inverso)
coefs_polyfit = np.polyfit(x, y, deg=2)
print(f"np.polyfit: {np.round(coefs_polyfit, 3)}")

# Usando la API moderna Polynomial
from numpy.polynomial import polynomial as P
coefs_fit = P.polyfit(x, y, deg=2)
print(f"Polynomial.polyfit: {np.round(coefs_fit, 3)}")

El ajuste por mínimos cuadrados con lstsq es la base de la regresión lineal y de muchas técnicas de modelado. La ventaja de formularlo como un sistema matricial es que se extiende de forma natural a múltiples variables predictoras y a términos no lineales.