1. ¿Qué es la matriz Jacobiana?

Para una función vectorial multivariable \(\mathbf{f}:\mathbb{R}^n \rightarrow \mathbb{R}^m\), la matriz Jacobiana \(J\) es la matriz de todas las derivadas parciales de primer orden:

\[J(\mathbf{x}) = \begin{bmatrix} \frac{\partial f_1}{\partial x_1} & \dots & \frac{\partial f_1}{\partial x_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial f_m}{\partial x_1} & \dots & \frac{\partial f_m}{\partial x_n} \end{bmatrix}\]

Esta matriz captura cómo cambian todas las salidas respecto a cada una de las entradas y es fundamental en:

• Optimización sin restricciones (Newton‑Raphson, métodos de quasi‑Newton).
• Control de sistemas no lineales y robótica.
• Back‑propagation en redes neuronales.
• Modelado de sensibilidad y análisis de incertidumbre.

2. Algoritmos clásicos para calcular el Jacobiano

3. Autodiferenciación: la solución moderna

Las librerías de autodiferenciación (AD) construyen el Jacobiano a partir del grafo computacional de la función, combinando precisión simbólica con velocidad de cálculo numérico.

• Forward‑mode AD: ideal cuando \(n \le m\).
• Reverse‑mode AD: la base del back‑propagation, óptimo para \(m \le n\) (p.ej., redes neuronales).

En Python, las bibliotecas más usadas son autograd, JAX, PyTorch y TensorFlow.

4. Implementaciones en Python

4.1. Con NumPy + diferencias finitas

import numpy as np
def jacobian_fd(fun, x, eps=1e-8):
    x = np.asarray(x, dtype=float)
    n = x.size
    f0 = np.asarray(fun(x))
    m = f0.size
    J = np.empty((m, n))
    for j in range(n):
        x_eps = x.copy()
        x_eps[j] += eps
        f1 = np.asarray(fun(x_eps))
        J[:, j] = (f1 - f0) / eps
    return J
# Ejemplo de uso
f = lambda v: np.array([v[0]**2 + v[1], np.sin(v[0]*v[1])])
print(jacobian_fd(f, np.array([1.0, 2.0])))

Esta rutina es útil para prototipos rápidos, pero recuerda calibrar eps según la escala de tus variables.

4.2. Con SymPy (cálculo simbólico)

import sympy as sp
x1, x2 = sp.symbols('x1 x2')
F = sp.Matrix([x1**2 + x2, sp.sin(x1*x2)])
J = F.jacobian([x1, x2])
print(J)
# Evaluar en un punto concreto
print(J.subs({x1: 1, x2: 2}).evalf())

Ideal cuando necesitas una expresión analítica para análisis de sensibilidad o generación de código C.

4.3. Con autograd (forward‑mode)

import autograd.numpy as anp
from autograd import jacobian
def f(x):
    return anp.array([x[0]**2 + x[1], anp.sin(x[0]*x[1])])
J = jacobian(f)
print(J(anp.array([1.0, 2.0])))

Sin necesidad de escribir diferencias finitas; la precisión es del orden de máquina.

4.4. Con PyTorch (reverse‑mode, ideal para redes)

import torch
def f(x):
    return torch.stack([x[0]**2 + x[1], torch.sin(x[0]*x[1])])
x = torch.tensor([1.0, 2.0], requires_grad=True)
# Forward pass
y = f(x)
# Jacobiano: grad(y_i) w.r.t x para cada i
J = []
for i in range(y.shape[0]):
    grad_i = torch.autograd.grad(y[i], x, retain_graph=True)[0]
    J.append(grad_i)
J = torch.stack(J)
print(J)

Esta técnica se usa internamente en el entrenamiento de modelos profundos.

5. Comparativa de rendimiento (Python 3.11, CPU Intel i7‑12700)

6. Buenas prácticas y troubleshooting

6.1. Elección del método

• Si la función está escrita en Python y ya usa torch o jax, aprovecha la autodiferenciación integrada.
• Para modelos de simulación en C/C++ que no pueden modificarse, prefiere diferencias finitas con scipy.optimize.approx_fprime.

6.2. Problemas comunes

• Valores NaN o Inf en el Jacobiano: suele deberse a puntos donde la función no es diferenciable o a overflow en operaciones trigonométricas. Use np.clip o regularice la función.
• Jacobiano escaso: Cuando la mayoría de derivadas son cero, convierta la matriz a formato scipy.sparse.csr_matrix para ahorrar memoria.
• Inestabilidad numérica en diferencias finitas: ajuste eps con la regla de O'Leary: \(\epsilon = \sqrt{\text{machine\_eps}} \times (1+|x|)\).

6.3. Seguridad y aislamiento

Si el cálculo del Jacobiano se ejecuta con datos provenientes de usuarios externos (p. ej., en una API de optimización), valide que los inputs sean numéricos y de tamaño razonable para evitar ataques de denegación de servicio mediante vectores de dimensión explosiva.

7. Casos de uso reales

• Robótica: cálculo de la Jacobiana del manipulador para controlar velocidad y fuerza en el espacio de trabajo.
• Optimización de hiperparámetros: métodos basados en Newton requieren el Hessiano, que es la derivada del Jacobiano.
• Modelos físicos diferenciables: simulaciones de fluidos donde la fuerza depende de parámetros aprendidos mediante AD.