Gauss-Newton en Python. Mínimos cuadrados no-lineales n-dimensionales.

El problema de optimización no lineal mediante mínimos cuadrados (no sujeto a restricciones) se enuncia formalmente como:
\displaystyle \underset{x}{minimizar}\, f(x)=\underset{i=1}{\overset{m}{\sum}}f_{i}(x)^{2}\equiv\frac{1}{2}F(x)^{T}F(x)
donde la función objetivo está definida en términos de funciones auxiliares {\left\{ f_{i}\right\} } y {F(x)} es una función vectorizada:
\displaystyle F(x)=(f_{1}(x)\: f_{2}(x)\:\cdots\: f_{m}(x))^{T}
En el contexto del ajuste de datos, las funciones auxiliares {\left\{ f_{i}\right\} } no son funciones arbitrarias, sino que corresponden los residuales del ajuste.
Las componentes del gradiente {\nabla f(x)} pueden calcularse como sigue:
\displaystyle \nabla f(x)=\nabla F(x)F(x)
De la misma forma el valor de la matriz Hessiana:
\displaystyle \nabla^{2}f(x)=\nabla F(x)\nabla F(x)^{T}+\underset{i=1}{\overset{m}{\sum}}f_{i}(x)\nabla^{2}f_{i}(x)
La matriz Hessiana es en muchas ocasiones complicada, cuando no imposible, de obtener. Tendremos que realizar {n^{2}} derivadas parciales dobles, por lo que para un problema como el que sigue de 5 variables tendremos que calcular 25 de estas operaciones. En cualquier caso su cálculo es pesado tanto si lo realizamos analíticamente como si aproximamos mediante diferencias finitas, aquí la potencia del método:
Llamemos {x_{*}} a la solución del problema, que causaría un ajuste perfecto de los datos, o lo que es lo mismo:
\displaystyle F(x)\approx0\: para\, x\approx x_{*}
Lo que nos lleva a:
\displaystyle \nabla^{2}f(x)=\nabla F(x)\nabla F(x)^{T}+\underset{i=1}{\overset{m}{\sum}}f_{i}(x)\nabla^{2}f_{i}(x)\approx\nabla^{2}f(x)=\nabla F(x)\nabla F(x)^{T}
\displaystyle
Por lo tanto este método sólo requiere el cálculo de las primeras derivadas para el cálculo del gradiente. Ahora sólo tendremos que resolver el sistema de ecuaciones alegraicas:
\displaystyle \nabla^{2}f(x)p=-\nabla f(x)
Implementación del algoritmo.
Mi implementación es un borrador para resolver un problema de optimización que consiste en lo siguiente, sólo para poneros en contexto. debo encontrar los valores para las constantes {H,\sigma,B,x_{o},y_{o}} que minimizen la distancia de la superficie definida por una campana doble de Gauss:
\displaystyle G=He^{-((x-x_{0})^{2}+(y-y_{0})^{2})\frac{1}{\sigma^{2}}}+B
respecto de los datos de brillo de una estrella, los cuales obtengo de una imagen directamente (podremos considerarlos a todos los efectos datos del problema).
Quien juegue un poco con numpy en el ámbito de la optimización puede venirle a la cabeza algunas funciones que proporciona, como la de mínimos cuadrados. No hay ningún problema con ellas, aunque no se si pueden trabajar con arrays n-dimensionales o mallas, pero el caso es que necesitaba tener control sobre el algoritmo y poder cambiar algunas partes.
Los pasos para resolver el algoritmo son los siguientes:
  • Hacer una estimación inicial de la semilla a la que llamaré {\beta}.
  • Calcular el gradiente, nabla, de la función G definida anteoriormente. Este calculo se hará analíniticamente porque siempre necesito aplicar la misma función.
  • Calcular la matriz Hessiana llamada “Hessian”.
  • Resolver el sistema de ecuaciones algebraicas para obtener unos nuevos pesos.
  • Sumar el vector de pesos a beta y realizar una nueva iteración hasta satisfacer el criterio de parada.
Vamos ahora con el código:
def Gaussian_distrib_3d(beta,x,y):
H,B,sigma,x0,y0 = beta
exponente = -( (x-x0)**2 +(y-y0)**2) / (2*sigma*sigma)
return H*np.e**(exponente) + B
Ahora vamos con el cuerpo del algoritmo:
gradientF = np.zeros([Nvariables,Ndata])

for step in range(0,N):
 H,B,sigma,x0,y0 = beta
 # Calculo de algunos parametros que se repiten
 p1 = (meshx-x0)**2 +(meshy-y0)**2
 p2 = 2*sigma*sigma
 p3 = sigma**3
 p4 = sigma*sigma
 exponente = -p1/p2

 Fx = np.ravel( H*np.e**(exponente) + B -z )
 Fx = Fx.transpose()

 # Calculo del gradiente
 gradientF[0,:] = np.ravel(np.e**(exponente))
 gradientF[1,:] = 1
 gradientF[2,:] = np.ravel((p1/p3)*H*np.e**(exponente))
 gradientF[3,:] = np.ravel(((meshx-x0)/p4) * H * np.e**(exponente))
 gradientF[4,:] = np.ravel(((meshy-y0)/p4) * H * np.e**(exponente))

 gradientFTransp = gradientF.transpose()

 nabla = np.dot(gradientF,Fx).reshape(Nvariables,1)

 #Calculamos la aproximacion de la matriz Hessiana
 Hessian = np.dot(gradientF,gradientFTransp)

 #Resolvemos el sistema algebraico
 p = np.linalg.solve(Hessian, -nabla)
 beta += p
El algoritmo se inicia tomando unos valores iniciales para la semilla, beta en nuestro caso:
beta = np.array([100,6,1.2,12,10.5]).reshape(Nvariables,1)
Por último los valores de x e y, puesto que se trata de una función ‘tridimensional’ (no en términos matemáticos) los generaremos mediante la función meshgrid()
z = data #datos que buscamos aproximar.
x = np.arange(24)
y = np.arange(24)
meshx,meshy = np.meshgrid(x,y)

Este algoritmo no es tan muy a la aproximación inicial, pero sí al conjunto de datos que buscamos ajustar. Si los residuales son muy altos o hay datos fuera de rango, ‘outers’, el algoritmo será inestable y seguramente lleve a una matriz singular que no podremos resolver directamente y deberemos resolverla por mínimos cuadrados. Esto último no será en la mayoría de los casos una opción viable.

El algoritmo toma pocas iteraciones para ajustar con una tolerancia muy baja. Un resultado para un subcaso bidimensional del problema anterior que empleé para validar el algoritmo:

fit

En el caso tridimensional el algoritmo funciona igualmente bien;

3adjust

Finalmente el mapa de errores, como se puede ser el ajuste es bastante bueno.

error_map

Conclusiones.

No estoy seguro de si las herramientas de las que dispone numpy pueden manejar los datos de la forma en que mi problema los dispone, la función numpy.optimize.lsqr() es versátil y seguramente pueda manejar arrays generados por la función meshgrid() si los pasamos propiamente aplanados (mediante ravel() o reshape() ).
En cualquier caso es un algoritmo fácil de implementar para problemas de optimización complejos y de varias variables en los que tenemos más o menos claro cuál es el tipo de función que va a ajustar bien con los datos.
Puesto que el factor sensible del algoritmo son los datos fuera de rango es recomendable preprocesar los mismos mediante alguna herramienta estadística (en breve pondré algo relacionado).
Esta entrada fue publicada en matemáticas, python y etiquetada , , , , , , , . Guarda el enlace permanente.

3 respuestas a Gauss-Newton en Python. Mínimos cuadrados no-lineales n-dimensionales.

  1. Pingback: Uso de numpy.optimize.leastsq() o cómo ahorrar trabajo. | Por las barbas de Kutta

  2. Juanlu001 dijo:

    Te tengo que confesar que de ajuste paramétrico ahora mismo no tengo ni idea y que tu entrada me ha producido una curiosidad enorme😀 ¡Estupenda entrada!

    En el código al principio me preguntaba por qué usabas `ravel` todo el rato, pero ya veo que es para generalizar el algoritmo a 2 dimensiones.

    Poco tengo que añadir, lo único ¿podrías publicar el código entero? A primera vista no aprecio muy bien cómo has integrado los bloques.

    ¡Un saludo!

    Juanlu

    • Samuel dijo:

      ¡Gracias! Es un tema muy interesante!

      No quise poner todo el código por no aburrir con muchas líneas y porque pensé que al no tratarse de un caso generalizado quizás liaría más que aclararía. Lo que voy a hacer es poner un enlace al notebook (si averiguo cómo) para que se pueda ver el proceso completo. Gracias por el feedback!

      Edit: la finalidad de ravel() ya la has visto, es un poco chapuza porque visto con más calma puedo hacer el ravel sobre meshx y meshy y generar directamente el gradiente en la forma ‘traspuesta’, si traspongo este tendría el original y así me ahorraría todo ese cálculo dentro del bucle… es lo malo de escribir código sobre la marcha sin pensar demasiado!

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s