Aproximación de raíces, puntos fijos y optimización de funciones

In [1]:
f = @(x)(1 + (x.^3 - 4*x) + log(1+x.^2));

xs = linspace(-3.5,2.5,100);
ys = f(xs);
plot(xs,ys)

Resumen: raíces y puntos fijos en una variable

  • Si f(a) y f(b) tienen signos opuestos $\Rightarrow$ f tiene una raíz en el intervalo [a,b]
    • Método de bisección, que produce cotas del error cometido.
  • Si $g:[a,b]\rightarrow [a,b]$ $\Rightarrow$ g tiene un punto fijo en el intervalo [a,b]
  • Si $g:[a,b]\rightarrow [a,b]$ y $|g'(x)|\leq k < 1$, $\Rightarrow$
    • $g$ tiene un único punto fijo $x^*$ en el intervalo [a,b]
    • la sucesión $x_{n+1}=g(x_n)$ converge a $x^*$, para cualquier punto inicial $x_0$.
    • $|x_n - x^*| < \frac{k^n}{1-k}|x_1-x_0|$

Resumen: métodos de Newton y Secante en una variable

  • El método de Newton busca raíces de f(x), iterando $$x_{n+1}=x_{n} - \frac{f(x_n)}{f'(x_{n})}$$
    • Es una iteración de punto fijo, con $g(x) = x - f(x)/f'(x)$.
    • No es fácil saber si convergerá... pero cuando lo hace es de orden 2
    • Hay que evaluar la derivada
  • El método de la secante es similar al de Newton, salvo porque: $$x_{n+1}=x_{n} - \frac{f(x_n)\left(x_{n}-x_{n-1}\right)}{f(x_{n})-f(x_{n-1})}$$
    • No hay que evaluar la derivada
    • El orden de convergencia es algo menor

Curiosidad: Mostramos en qué raíz termina el método de Newton (cuando termina), dependiendo del punto de inicio.

cuencas

Métodos híbridos

Combinando métodos como los anteriores, se obtienen métodos prácticos como el de Brent.

Ejercicio:

In [ ]:
In [ ]:

Ejercicio:

En la documentación oficial de matlab para fzero, averigüa:

  • Cómo obtener el número de iteraciones del algoritmo y el número de evaluaciones de la función.
  • Cómo modificar la tolerancia en la x.
  • Cómo modificar la tolerancia en la y.
In [ ]:
In [ ]:

Raíces y puntos fijos en varias variables

  • Si f(a) y f(b) tienen signos opuestos $\Rightarrow$ f tiene una raíz en el intervalo [a,b]

Tiene un equivalente teórico en varias variables, pero no tiene aplicación práctica general.

  • Si $g:[a,b]\rightarrow [a,b]$ $\Rightarrow$ g tiene un punto fijo en el intervalo [a,b]

En varias variables, teorema del punto fijo de Brouwer:

  • Toda aplicación continua $g:\mathbf{B}\rightarrow\mathbf{B}$, donde $\mathbf{B}\subset\mathbb{R}^n$ es un conjunto "sin agujeros" (por ejemplo, si es convexo), tiene al menos un punto fijo.
  • Si $g:[a,b]\rightarrow [a,b]$ y $|g'(x)|\leq k < 1$, $\Rightarrow$
    • $g$ tiene un único punto fijo $x^*$ en [a,b]
    • $x_{n+1}=g(x_n)$ converge a $x^*$, para cualquier punto inicial $x_0$.
    • $|x_n - x^*| < \frac{k^n}{1-k}|x_1-x_0|$

Cierto, casi sin cambios, en varias variables:

  • Si $g:\mathbf{B}\rightarrow \mathbf{B}$ y $\Vert D g(x)\Vert\leq k < 1$, $\Rightarrow$
    • $g$ tiene un único punto fijo $x^*$ en $\mathbf{B}$
    • la sucesión $x_{n+1}=g(x_n)$ converge a $x^*$, para cualquier punto inicial $x_0$.
    • $|x_n - x^*| < \frac{k^n}{1-k}|x_1-x_0|$

Si $g:\mathbb{R}^n\rightarrow\mathbb{R}^n$, tenemos, por el teorema fundamental del cálculo, para $y=x+h$:

$$ g(y)-g(x)=\left(\int _{0}^{1}Dg(x+th)\cdot h\,dt\right)$$

luego si $\Vert D g(x) \Vert \leq k$ para todo $x$, tenemos que

$$ \begin{split} \Vert g(y) - g(x)\Vert &=\left\Vert\int _{0}^{1}Dg(x+th)\cdot h\,dt\right\Vert\\ &\leq \int _{0}^{1} \Vert Dg(x+th)\cdot h\Vert\,dt\\ &\leq k\int _{0}^{1} \Vert h\Vert\,dt \leq k \Vert y-x\Vert \end{split} $$

Por tanto, si $\Vert D g(x)\Vert\leq k < 1$, aplicar g reduce distancias.

Nota: No hemos hablado aún de las normas matriciales, $\Vert D g(x) \Vert$, pero veremos más adelante que para una norma adecuada, $\Vert D g(x) \Vert \leq k\Rightarrow \Vert D g(x)\cdot v \Vert \leq k\Vert v\Vert$.

Una aplicación tal que $\Vert g(y)-g(x)\Vert \leq k\Vert y- x\Vert$, para $k<1$, se llama contractiva, y la imagen de un conjunto se reduce hasta que converge a un punto.

contractive mapping

Definimos la sucesión $x_{n+1} = g(x_n)$, comenzando con un punto cualquiera $x_0$:

$$\Vert x_{n+1} - x_n\Vert = \Vert g(x_n) - g(x_{n-1})\Vert \leq k \Vert x_n-x_{n-1}\Vert $$

y del mismo modo que en una variable:

$$\Vert x_n - x^*\Vert < k^n\Vert x_0-x^\ast\Vert$$

o tb:

$$\Vert x_n - x^*\Vert < \frac{k^n}{1-k}\Vert x_1-x_0\Vert$$

Relación entre optimización y minimización

Los problemas de encontrar raíces y de minimizar (o maximizar) funciones son parecidos: $$f:\mathbb{R}^n\rightarrow\mathbb{R}^n$$ $$h:\mathbb{R}^n\rightarrow\mathbb{R}$$

  • Los mínimos de $h$ son puntos críticos $\Rightarrow$ son soluciones del sistema $\nabla h = \mathbf{0}$.
  • Si $h$ es convexa, el único mínimo de $h$ es la única raíz de $\nabla h = \mathbf{0}$.
  • Una función puede tener muchas raíces en un intervalo, y tb muchos mínimos y máximos locales.
  • $f(\mathbf{x_0})=\mathbf{0}$ $\Rightarrow$ $\Vert f(\mathbf{x})\Vert^2$ alcanza un mínimo absoluto en $\mathbf{x_0}$.
  • Hay un análogo directo al método de bisección para minimizar funciones de una variable: Golden-ratio search

En funciones convexas, minimizar f y encontrar raíces de f' es equivalente (incluso "si f no es diferenciable").

Sin embargo, hay una diferencia importante:

  • Todas las ráices de una función son equivalentes a priori: nos basta con una, o bien las queremos todas.
  • Cuando buscamos el mínimo de una función, nos interesa sobre todo el mínimo absoluto. A menudo es preferible un punto x con valor pequeño de f a un mínimo local con mayor valor de f (f(y)>f(x)).

Teorema de Weierstrass

Una función continua siempre alcanza su valor máximo y un mínimo en un conjunto "cerrado" y acotado.

También conocido como Teorema de los valores extremos

Problemas:

  • Una función definida en un conjunto no acotado puede no tener un mínimo ni un mínimo.
  • El máximo o el mínimo puede estar en el interior, o en la frontera del conjunto...
  • Una función puede tener muchos máximos y mínimos locales

Mínimos locales y globales

Si una función tiene muchos mínimos locales, un algoritmo de búsqueda de mínimos puede "quedarse atascado" en un mínimo local, y no alcanzar el mínimo absoluto...

In [8]:
f = @(x)((x - 2) .* x .* (x + 2) .* (x + 3));

xs = linspace(-3.5,2.5,100);
ys = f(xs);

plot(xs,ys)
In [5]:
res1 = fminbnd(f, -3, 0);
res2 = fminbnd(f,  0, 2);
minima = [res1, res2];
f = @(x)((x - 2) .* x .* (x + 2) .* (x + 3));

xs = linspace(-3.5,2.5,100);
ys = f(xs);

hold on
plot(xs,ys)
plot(minima, f(minima), 'og')
hold off

Ejercicio

  • Investiga la diferencia entre las funciones fminbnd y fminsearch en la documentación oficial de matlab.
  • Investiga cómo obtener información adicional:
    • El valor de f en el punto
    • El número de evaluaciones de f que han sido necesarias
  • Investiga cómo fijar la tolerancia.
In [12]:
fminsearch(f, -4)

ans = -2.5742
In [13]:
fminsearch(f, -2)

ans =  1.2538
In [27]:
def f(x):
    return (x - 2) * x * (x + 2)* (x+3)*(1+(x-1)**2)

res1 = fminsearch(f, -4);
print(res1.x, res1.fun)

res2 = fminbnd(f, -3, -1);
print(res2.x, res2.fun)

1.492848664616281 -14.766957914752735
-2.629300965570816 -40.240314518277174
In [32]:
xs = np.linspace(-3.5,2.5,100)
ys = f(xs)
fig, ax = plt.subplots(figsize=(9,5))
ax.plot(xs,ys)
ax.plot([res1.x, res2.x], [res1.fun, res2.fun], 'go')
ax.set_ybound(-50,50)
plt.show()

Método de Newton en varias variables

  • El método de Newton busca raíces de f(x), iterando $$x_{n+1}=x_{n} - \frac{f(x_n)}{f'(x_{n})}$$

  • El método de Newton en varias variables usa esta iteración: $$x_{n+1}=x_{n} - \left(\nabla f(x_n)\right)^{-1}f(x_n)$$

    • Es una iteración de punto fijo, con $g(x) = x - \left(\nabla f(x_n)\right)^{-1}f(x)$.
    • No es fácil saber si convergerá... pero cuando lo hace es de orden 2
    • Hay que evaluar la matriz jacobiana, pero no es necesario invertirla, basta con resolver el sistema $\left(\nabla f(x_n)\right)y_n = f(x_n)$

Más métodos numéricos para encontrar raíces

Para encontrar raíces de funciones $f:\mathbb{R}^n\rightarrow\mathbb{R}^n$ hay muchos métodos (matlab):

scipy.optimize.root (scipy.optimize)

    method : str, optional
    Type of solver. Should be one of

            ‘hybr’ ‘lm’ ‘broyden1’ ‘broyden2’ ‘anderson’ ‘linearmixing’ ‘diagbroyden’ ‘excitingmixing’ ‘krylov’ ‘df-sane’

Para encontrar raíces de funciones $f:\mathbb{R}^n\rightarrow\mathbb{R}^n$ hay muchos métodos:

  • Ninguno ofrece garantías completas de encontrar la raíz.
  • Algunos de hecho intentan minimizar $\Vert f(x)\Vert^2$, y reportan un resultado incluso cuando $f(x)=0$ no tiene solución.
  • Los de tipo "Broyden", son similares al de Newton, pero intentan evitar calcular el gradiente (en el espíritu del método de la secante).
  • Ejercicio: Averiguad qué método usa matlab (por defecto) para resolver sistemas de ecuaciones no lineales.

Optimización de funciones de variables variables

  • Sabemos que es imposible encontrar el valor óptimo de forma garantizada, excepto si f es convexa.
  • La optimización puede terminar en un mínimo o máximo local.

Algunas ideas importantes:

Descenso del gradiente y sus variantes:

Para acercarte a un mínimo de $f$, avanza un poco en la dirección $-\nabla f(x)$

  • Los algoritmos de esta familia tardan poco en alcanzar un mínimo local.
  • El mínimo local puede tener un valor muy superior al mínimo global.

Simmulated annealing

(el objetivo es maximizar la función)

  • Una temperatura alta acepta soluciones inferiores al "mejor valor encontrado hasta el momento".
  • Una temperatura baja no acepta soluciones inferiores al "mejor valor encontrado hasta el momento".
  • Comienza con una temperatura alta, y la baja progresivamente...

Hill_Climbing_with_Simulated_Annealing

Basin hopping

  • Un salto aleatorio
  • Después sigue el gradiente hasta un mínimo local
  • Acepta o rechaza el cambio en función de la temperatura

Genetic algorithms

Combina varias ideas de inspiración biológica:

  • Almacena varios puntos $x_i$, en lugar de sólo uno
  • Mutaciones
  • Recombinaciones
  • Selección natural

... y muchos más

'Nelder-Mead' 'Powell' 'CG' 'BFGS' 'Newton-CG' 'L-BFGS-B' 'TNC' 'COBYLA' 'SLSQP' 'trust-constr' 'dogleg' 'trust-ncg' 'trust-krylov' 'trust-exact'

Ejercicios

  • Probad varios métodos para encontrar raíces de funciones en varias variables.
  • Probad varios métodos para minimizar funciones en varias variables.

    ¡No es necesario que lo programéis!

  • Responded a estas preguntas:
    • ¿por qué las librerías de software ofrecen métodos de minimización pero no de maximización?
    • ¿qué métodos implementan las funciones que has elegido? ¿necesitas conocer el gradiente de la función?
    • ¿hay alguna circunstancia en la que puedas asegurar que los métodos que has elegido funcionarán?
In [ ]:

Ejercicio

Queremos determinar el flujo de agua a través de una tubería.

Encuentra el factor de fricción de Darcy que resuelve la ecuación de Colebrook - White

$${\displaystyle {\frac {1}{\sqrt {f}}}=-2\log \left({\frac {\varepsilon }{3.7D_{\mathrm {h} }}}+{\frac {2.51}{\mathrm {Re} {\sqrt {f}}}}\right)} $$

con los siguientes datos:

  • Diámetro $D_h$: 8e-2
  • Número de Reynolds: 2e5
  • Rugosidad $\varepsilon$: 2e-4
In [ ]:
In [ ]:
In [ ]:

Ejercicio

Encuentra una solución, si la hay, del sistema de ecuaciones: $$ \begin{array}{rr} x^2 - y^2 &= 2\\ xy&=2 \end{array} $$ Prueba al menos un par de métodos distintos. ¿Encuentras siempre la misma solución?

In [ ]:
In [ ]:
In [ ]:

Ejercicio

Busca un mínimo global de la función de Rosenbrock de dos variables:

$$f(x,y)=(a-x)^{2}+b(y-x^{2})^{2}$$

para a=1, b=100, usando al menos dos métodos distintos.

In [2]:
In [ ]:
In [ ]:

Ejercicio

Consideramos la aplicación $F(x,y) = \left(\frac{y}{2}, \frac{-x}{2} + 10 \right)$. El profesor te garantiza que es contractiva con constante $k=\frac{1}{2}$ (es la composición de una rotación, una homotecia de factor 1/2 y un desplazamiento).

  • Calcula las iteraciones sucesivas de esta función partiendo de distintos puntos iniciales.
  • Dibuja la evolución del error cometido como función del número de iteraciones.
  • Compara las estimaciones del error con las vistas más arriba.
In [ ]:
In [ ]:
In [ ]:

Ejercicio

Repite el ejercicio anterior para $$F(x,y) = \left(\cos\left(\frac{y}{2}\right)+5, \frac{x+y}{3} + 1 \right)$$

  • Calcula las iteraciones sucesivas de esta función partiendo de distintos puntos iniciales.
  • Dibuja la evolución del error cometido como función del número de iteraciones.
  • Compara las estimaciones del error con las vistas más arriba.
  • ¿Crees que es contractiva?
In [ ]:
In [ ]:
In [ ]: