H ay un momento, durante el grado, en el que te das cuenta de que casi todo lo que vas a hacer como ingeniero es resolver la misma ecuación una y otra vez. Calor, electromagnetismo, fluidos, sonido — al final son variaciones de una EDP de segundo orden con condiciones de contorno horribles. Y para todas ellas, alguien ya programó un solver mejor que cualquier cosa que vayas a escribir tú esa tarde.
Pero hay un agujero. Un caso donde los solvers clásicos tropiezan: cuando no conoces los parámetros del sistema y solo tienes datos dispersos.
Ahí entran los Physics-Informed Neural Networks. Este post va de eso, pero contado por las cosas que me funcionaron y las que no.
Lo que dejó de cuadrarme con FEM
En clase de vibraciones nos contaron FEM como el martillo universal. Mallar, ensamblar, resolver. Bonito, riguroso, con teoría detrás. Funciona.
Luego empecé a hacer cosas por mi cuenta y aparecieron las pequeñas miserias. Mallar una geometría irregular en 3D — un conducto, un altavoz con su recinto, un detector ultrasónico — lleva más tiempo que la simulación. A frecuencias altas necesitas entre 6 y 10 elementos por longitud de onda, lo que en aire a 15 kHz son unos 3 mm por elemento. Una sala normal te pide millones.
FDTD tiene un atractivo distinto, casi infantil: una rejilla, una receta, un bucle. Lo programé un par de veces en cosas de juguete (en Python, ineficiente) y la primera vez que metí 3D el portátil me dijo claramente que no. La condición CFL te aprieta los pasos de tiempo y la memoria te aprieta todo lo demás.
Hasta aquí, problemas conocidos y abordables. La pesadilla real es otra: el problema inverso.
Tienes mediciones de presión acústica en cinco micrófonos. ¿Cuál era la fuente? ¿Cómo varía la velocidad del sonido c(x) dentro del medio (que no es homogéneo, claro, porque la vida no es homogénea)? Cada evaluación pide resolver la PDE entera. Para derivar respecto a c necesitas el método adjoint — otra PDE, hacia atrás en el tiempo, acoplada. Funciona, pero es lento, frágil, y los resultados oscilan a no ser que ajustes la regularización con cariño.
La primera vez que vi a unos investigadores reportar el mismo problema resuelto con un PINN en cinco páginas y un script de PyTorch, no me lo creí. Luego lo probé.
La idea, sin envoltorio
Una red neuronal p_θ(x, t) es un aproximador universal con derivadas exactas respecto a sus entradas (gracias a autodiff). La idea de Raissi, Perdikaris y Karniadakis (2019) cabe en una línea:
Si conoces la PDE que el campo debe cumplir, métela dentro de la función de pérdida.
No discretizas el dominio. Muestreas puntos (xᵢ, tᵢ) aleatorios (o cuasi-aleatorios, que va mejor), evalúas la red ahí, calculas las derivadas con torch.autograd.grad, y sumas el cuadrado del residuo de la PDE como una pérdida más. El optimizador empuja los pesos hacia el cumplimiento simultáneo de datos, física y contornos.
La primera vez que lo entendí pensé que era un truco barato. Al fin y al cabo, una red neuronal puede aprender cualquier cosa si le das datos suficientes. Lo brillante es que funciona también cuando los datos son escasos o no hay. La física hace el trabajo pesado y los pocos datos disponibles afinan la solución.
La ecuación de onda, sin maquillaje
Residuo:
Y la pérdida de PDE como la media de r_θ² sobre N puntos de colocación:
Sumas contornos e iniciales, y a entrenar. Fácil de escribir. Lo difícil viene en el detalle, y es donde más tiempo perdí.
El sesgo espectral, o por qué mi primera onda parecía gelatina
Mi primer PINN para ondas 1D fue un desastre instructivo. La red — MLP de cuatro capas con tanh — aprendió una "envolvente" suave del campo y se negó tercamente a aprender las oscilaciones. Salía algo que parecía la presión promediada en el tiempo. Pasé dos tardes convencido de que tenía un bug en el código.
No era un bug. Es un fenómeno conocido: sesgo espectral. Las redes feedforward estándar aprenden frecuencias bajas mucho antes que altas. Para acústica, donde la información vive en frecuencias medias-altas, es letal.
Dos arreglos que funcionan:
Fourier features (Tancik et al., 2020). Proyectas la entrada con sin(2πBx), cos(2πBx) con B muestreada de una gaussiana. Cuesta tres líneas de código y cambia todo.
SIREN (Sitzmann et al., 2020). Activaciones senoidales en lugar de tanh. Más radical, más sensible a la inicialización.
Yo siempre arranco con Fourier features. Más fácil de domar.
Los λᵢ no son inocentes
Pones λ_PDE = 1, λ_BC = 1, λ_IC = 1 y rezas. Mal. Los gradientes de cada término viven en escalas distintas y casi siempre alguno aplasta a los demás. La red termina ajustando solo lo que más pesa.
Solución casera: ponderar a ojo. Yo arranco con λ_BC = λ_IC = 10 y λ_PDE = 1, y muchas veces basta. Solución elegante: pesos auto-adaptativos (Wang, Teng y Perdikaris, 2021), que se reescalan dinámicamente según los gradientes de cada término. Cuando tengas un problema serio, te ahorra horas.
La onda que se moría en el origen
Otra que me comí entera. Si entrenas con puntos de colocación repartidos uniformemente en el espacio-tiempo, el optimizador encuentra cómodo un mínimo local en el que la solución muere cerca del origen y deja el resto del dominio plano. Es la red haciendo trampa: cumple aproximadamente todo, no destaca en nada.
La solución es causal training (Wang, Sankaran y Perdikaris, 2022): ponderar la pérdida en t = tₖ por un factor que solo crece cuando la pérdida en t < tₖ ya está controlada. La red aprende en orden temporal, igual que hace la física. Una vez lo metí en mi código, los pulsos empezaron a viajar de verdad.
El sábado que me lo creí: problemas inversos
Aquí va la parte que llevó este post a existir.
Un sábado por la tarde, sobremesa, me dio por probar el modo inverso. Monté el ejemplo más estúpido que se me ocurrió: ecuación de onda 1D, velocidad real c = 1.5, condición inicial gaussiana, contornos de Dirichlet. Generé datos sintéticos muestreando la solución analítica en 30 puntos aleatorios de (x, t). Eso era todo lo que la red podía ver.
El cambio respecto al modo directo cabía en dos líneas:
c = torch.nn.Parameter(torch.tensor(0.5)) # ¡valor inicial deliberadamente malo!
optimizer = torch.optim.Adam(list(net.parameters()) + [c], lr=1e-3)
Sin penalización extra, sin regularizador, sin nada. La pérdida sumaba: residuo de PDE + condiciones iniciales + ajuste a los 30 puntos medidos. Le di a run y me fui a tomar un café.
Cuando volví, c estaba en 1.498. Tres milésimas de la verdad. Con 30 puntos. Sin haber resuelto explícitamente ninguna PDE.
Sé que suena a publicidad. Lo monté esperando que no funcionara y me sentí estúpido al ver que sí. La explicación es esta: el residuo de la ecuación de onda solo se anula simultáneamente sobre el campo correcto y con el c correcto. El optimizador no tiene otra ruta para bajar la pérdida que descubrir ambos. Es una de esas cosas que parecen magia y son matemática limpia.
Las aplicaciones reales de esta forma exacta del problema:
Caracterización no destructiva de materiales por ultrasonidos
Imagen ultrasónica cuantitativa (no solo modo B)
Full waveform inversion sísmica (misma ecuación, c del orden de km/s)
En audiología: recuperación de propiedades acústicas del conducto auditivo o de una prótesis a partir de mediciones in vivo
Lo último es lo que más me interesa, por razones que ya sabéis si seguís otros posts del blog.
El "hola mundo" en PyTorch
Versión limpia del experimento, sin la parte inversa para no liar. Dominio [0, 1] × [0, 1], c = 1, gaussiana inicial, contornos rígidos.
import torch
import torch.nn as nn
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
class PINN(nn.Module):
def __init__(self, hidden=64, depth=4, num_freq=16, sigma=5.0):
super().__init__()
self.B = torch.randn(2, num_freq, device=device) * sigma
layers = [nn.Linear(2 * num_freq, hidden), nn.Tanh()]
for _ in range(depth - 1):
layers += [nn.Linear(hidden, hidden), nn.Tanh()]
layers += [nn.Linear(hidden, 1)]
self.net = nn.Sequential(*layers)
def fourier(self, xt):
proj = 2 * torch.pi * xt @ self.B
return torch.cat([torch.sin(proj), torch.cos(proj)], dim=-1)
def forward(self, x, t):
xt = torch.stack([x, t], dim=-1)
return self.net(self.fourier(xt)).squeeze(-1)
def grad(y, x):
return torch.autograd.grad(y, x, grad_outputs=torch.ones_like(y),
create_graph=True)[0]
def pde_residual(model, x, t, c=1.0):
x.requires_grad_(True); t.requires_grad_(True)
p = model(x, t)
p_t = grad(p, t)
p_tt = grad(p_t, t)
p_x = grad(p, x)
p_xx = grad(p_x, x)
return p_tt - (c ** 2) * p_xx
model = PINN().to(device)
opt = torch.optim.Adam(model.parameters(), lr=1e-3)
N_pde, N_bc, N_ic = 5000, 200, 500
for step in range(20_000):
x = torch.rand(N_pde, device=device)
t = torch.rand(N_pde, device=device)
r = pde_residual(model, x, t)
loss_pde = (r ** 2).mean()
t_b = torch.rand(N_bc, device=device)
p0 = model(torch.zeros_like(t_b), t_b)
p1 = model(torch.ones_like(t_b), t_b)
loss_bc = (p0 ** 2).mean() + (p1 ** 2).mean()
x_i = torch.rand(N_ic, device=device); x_i.requires_grad_(True)
t_i = torch.zeros_like(x_i); t_i.requires_grad_(True)
p_i = model(x_i, t_i)
p_i_t = grad(p_i, t_i)
target = torch.exp(-100 * (x_i - 0.5) ** 2)
loss_ic = ((p_i - target) ** 2).mean() + (p_i_t ** 2).mean()
loss = loss_pde + 10 * loss_bc + 10 * loss_ic
opt.zero_grad(); loss.backward(); opt.step()
if step % 1000 == 0:
print(f"step {step:5d} | pde {loss_pde.item():.2e} "
f"bc {loss_bc.item():.2e} ic {loss_ic.item():.2e}")
Setenta líneas. En una GPU modesta, minutos. Si lo grafías tras entrenar, ves la gaussiana inicial dividirse en dos pulsos que rebotan en los contornos. Justo lo que predice d'Alembert. Recomiendo hacerlo: hay algo satisfactorio en ver una solución de PDE emerger de un optimizador estadístico.
Para convertirlo en problema inverso: las dos líneas de arriba y un término más en la pérdida con tus datos medidos. Cinco minutos de trabajo, una eternidad de implicaciones.
Antes de tirar tu FEM por la ventana
Sería deshonesto terminar aquí sin frenar el entusiasmo. Lo que los PINNs no hacen bien todavía:
Más lentos que FEM en problemas directos bien planteados con geometría sencilla. Si tu problema es ese, FEM le da una paliza. Sin discusión.
Convergencia mal entendida teóricamente. Cuando un PINN no converge, diagnosticarlo es más arte que ciencia. La sensación es de estar haciendo alquimia.
Frecuencias muy altas siguen siendo difíciles incluso con Fourier features.
Producción industrial todavía no. Hay esfuerzos serios (NVIDIA Modulus, entre otros), pero el ecosistema FEM/BEM lleva décadas de validación. Para certificar un audífono médico nadie va a aceptar un PINN sin años más de evidencia.
Dónde sí compensa volcarse: problemas inversos, medios heterogéneos complejos, fusión de datos experimentales con modelos físicos parciales (los famosos sistemas grey-box), y geometrías parametrizadas donde quieres soluciones para todo un rango de parámetros.
Dónde mirar a continuación
Variantes que vale la pena conocer aunque sea de oídas: XPINN descompone el dominio en subdominios con redes independientes y acopla en las interfaces (escala mejor); B-PINN son bayesianos y cuantifican incertidumbre, lo cual en aplicaciones médicas no es opcional.
Y un cambio de marco cercano: los operadores neuronales. En vez de aprender una solución, aprenden el operador que mapea condiciones iniciales y parámetros a soluciones. DeepONet (Lu et al., 2019) y FNO (Li et al., 2020) son los pilares. Para acústica esto es prometedor: una vez entrenado el operador, evaluar una nueva configuración cuesta lo que un forward pass.
Mi próximo experimento, cuando saque tiempo, es montar un FNO sobre HRTFs sintéticas y ver si consigo generar funciones de transferencia personalizadas a partir de geometría craneal aproximada. Si sale algo, lo cuento por aquí.
Referencias
Raissi, Perdikaris y Karniadakis (2019). Physics-informed neural networks. JCP.
Moseley, Markham y Nissen-Meyer (2020). Solving the wave equation with PINNs.
Wang, Sankaran y Perdikaris (2022). Respecting causality for training PINNs.
Tancik et al. (2020). Fourier features let networks learn high frequency functions.... NeurIPS.
Sitzmann et al. (2020). Implicit neural representations with periodic activation functions (SIREN). NeurIPS.
Rasht-Behesht et al. (2022). PINNs for wave propagation and full waveform inversion.
Si te animas a probarlo, dos consejos. Usa Fourier features desde el minuto cero. Y no te creas la primera curva de pérdida bonita que veas: yo me pasé media tarde celebrando una red que había aprendido a predecir cero en todas partes, porque la pérdida bajaba muy suave. La física no perdona los descuidos, ni siquiera cuando le pones una red delante.
Inicia sesión para votar este artículo.