Buenas a tod@s,
llevo unas semanas dándome de cabezazos y no consigo dar con la tecla. Les comento mi problema a ver si alguien me pudiera arrojar algo de luz.
La cuestión es que tengo un sistema de ecuaciones no lineales (3 para ser exactos), con tres incognitas. He implementado el método de Newton Raphson multivariable en python y funciona, pero para mi sorpresa dependiendo del valor de las funciones me he encontrado que la convergencia puede ser realmente lenta. Buscando encontré que con el algoritmo de
line search se puede optimizar el tamaño del paso para mejorar la convergencia. El problema viene aquí no sé si lo tengo bien implementado porque me salen cosas muy raras.
la función que tengo es la siguiente
$$\bar{\sigma}={\left(\frac{1}{2}\left({\left |{S_1 - S_2}\right |}^m+{\left |{S_2 - S_3}\right |}^m+{\left |{S_3 - S_1}\right |}^m\right)\right)}^\frac{1}{m}$$
Por no complicar el asunto ya que $$S_i$$ son valores complejos, aunque si lo necesitaseis añadiría todas la fórmulas, las 3 funciones no lineales las obtengo de tres casos particulares ya que de otro modo las derivadas serían monstruosas. Así pues $$f_1 = f(b, c)$$, $$f_2 = f(a, c)$$ y $$f_3=f(a, b)$$. Donde $$a, b, c$$ son las incógnitas a calcular por medio de Newton Raphson. Con estas tres ecuaciones construyo el vector $$B=\left[f_1, f_2, f_3\right]$$
El Jacobiano quedaría de la siguiente manera:
J=\begin{bmatrix}{\nabla{f_1(x)}}\\{\nabla{f_2(x)}}\\{\nabla{f_3(x)}}\end{bmatrix}=\begin{bmatrix}{\frac{{\partial f_1}}{{\partial a}}}&{\frac{{\partial f_1}}{{\partial b}}}&{\frac{{\partial f_1}}{{\partial c}}}\\{\frac{{\partial f_2}}{{\partial a}}}&{\frac{{\partial f_2}}{{\partial b}}}&{\frac{{\partial f_2}}{{\partial c}}}\\{\frac{{\partial f_3}}{{\partial a}}}&{\frac{{\partial f_3}}{{\partial b}}}&{\frac{{\partial f_3}}{{\partial c}}}\end{bmatrix}
Debido a los tres casos especiales se obtiene que $$\frac{{\partial f_1}}{{\partial a}}=\frac{{\partial f_2}}{{\partial b}}=\frac{{\partial f_3}}{{\partial c}}=0$$.
Para obtener el vector incremento $$P$$ para la iteración actual en el método de Newton-Raphson:
$$P=J^{-1}·(-B)$$
siendo los nuevos valores tal que:
$$\begin{pmatrix}{a1}\\{b1}\\{c1}\end{pmatrix}=P+\begin{pmatrix}{a0}\\{b0}\\{c0}\end{pmatrix}$$
Hasta aquí todo bien. Ahora empieza la tortura con el
line search. Si lo he entendido bien, consiste en estimar un tamaño inicial lo suficientemente grande en una dirección de búsque determinada para optimizar el paso de $$a_0 \rightarrow{a_1}$$.
La dirección de búsqueda deber ser tal que asegure el decaimiento de la función de búsqueda. Dicha función es del tipo:
$$\phi(\alpha) = f(X_k+\alpha{p_k})$$
donde $$X_k$$ es el vector de coeficientes actuales $$a0, b0, c0$$, $$\alpha$$ es el tamaño del paso y $$p_k$$ es la dirección de búsqueda que para asegurar el decaimiento $$p_k=-\nabla{f(X_k)}$$.
A la hora de implementar mi propia función de
line search en python tengo algunas dudas como: ¿Cómo defino el tamaño de paso inicial? La bibliografía establece un valor de paso de 1. Podría usarlo como valor por defecto, pero ¿hay alguna forma de calcular un primer valor?
El otro problema viene a la hora del criterio de convergencia de $$\alpha$$. Si me he enterado bien, se puede usar la condición de Armijo (1966) la cual dice que $$\alpha$$ se actualizará multiplicando por un valor $$\rho$$; tal que $$\alpha=\alpha \rho$$; hasta que se cumpla la condición:
$$f(X_k + \alpha{p_k})\geq{f(X_k) + c\alpha{\nabla{f(X_k)p_k}}}$$
donde $$c$$ y $$\rho$$ son parámetros del control de búsqueda y están el en rango $$(0, 1)$$.
El código que tengo en python es el siguiente:
def line_search(s11, s22, s33, J, B, P, a0, b0, c0, m, alpha_max=1.0):
# S11, S22, S33, m son parámetros para lad funciones f1, f2, f3
# J: Jacobiano
# B: Vector [f1, f2, f3] para a0, b0, c0
# P: Vector incremento [a1-a0, b1-b0, c1-c0]
alpha = 1.0 # initial step size
c = 0.5 # parameter for sufficient decrease condition
rho = 0.9 # parameter for backtracking
f0 = getPhiVector(s11, s22, s33, a0, b0, c0, m)
while(alpha > 1e-8):
a1 = a0 + alpha * P[0]
b1 = b0 + alpha * P[1]
c1 = c0 + alpha * P[2]
fnext = getPhiVector(s11, s22, s33, a1, b1, c1, m)
if(fnext[0] <= f0[0] - c * alpha * np.dot(J[0], J[0]) and
fnext[1] <= f0[1] - c * alpha * np.dot(J[1], J[1]) and
fnext[2] <= f0[2] - c * alpha * np.dot(J[2], J[2])):
return alpha
else:
alpha *= rho
return alpha_maxP.D. no sé si hay en LateX algún comando para insertar código.
El $$\alpha$$ devuelto es aplicado a la actualización de los coeficientes en el método de Newton Raphson tal que:
$$\begin{pmatrix}{a1}\\{b1}\\{c1}\end{pmatrix}=\alpha{P}+\begin{pmatrix}{a0}\\{b0}\\{c0}\end{pmatrix}$$
El caso es que no sé si estoy aplicando bien el
line search, la concición de convergencia de $$\alpha$$ o qué está pasando pues para nada mejora la convergencia.
Muchas gracias de antemano.
