Clase 03 (05/04)

Derivadas

Derivadas para una función en una variable

Método forward

Dados $f\in C^2$ y $h>0$, desarrollando Taylor centrado en $x$ tenemos que: $$f(x+h)=f(x) + f'(x)h + \color{red}{\dfrac{f''(\xi)h^2}{2}} \qquad \xi\in(x,x+h)$$ $$\begin{array}{rll}\Rightarrow f'(x) & = & \dfrac{f(x+h)-f(x)}{h} - \color{red}{\dfrac{f''(\xi)h}{2}} \\ & = & \dfrac{f(x+h)-f(x)}{h} - \color{red}{O(h)} \end{array}$$

Luego, $f'(x)$ puede aproximarse como: $$f'(x) \approx \dfrac{f(x+h)-f(x)}{h}$$

Método Backward

Análogo al método forward. Dados $f\in C^2$ y $h>0$, desarrollando Taylor centrado en $x$ tenemos que: $$f(x-h)=f(x) - f'(x)h + \color{red}{\dfrac{f''(\xi)h^2}{2}} \qquad \xi\in(x-h,x)$$ $$\begin{array}{rll}\Rightarrow f'(x) & = & \dfrac{f(x)-f(x-h)}{h} + \color{red}{\dfrac{f''(\xi)h}{2}} \\ & = & \dfrac{f(x)-f(x-h)}{h} + \color{red}{O(h)} \end{array}$$

Luego, $f'(x)$ puede aproximarse como: $$f'(x) \approx \dfrac{f(x)-f(x-h)}{h}$$

Método de centradas

Dados $f\in C^3$ y $h>0$, desarrollando Taylor centrado en $x$ tenemos que: $$f(x+h)=f(x) + f'(x)h + \dfrac{f''(x)h^2}{2} + \color{red}{\dfrac{f'''(\xi)h^3}{6}} \qquad \xi\in(x,x+h)$$ $$f(x-h)=f(x) - f'(x)h + \dfrac{f''(x)h^2}{2} - \color{red}{\dfrac{f'''(\eta)h^3}{6}} \qquad \eta\in(x-h,x)$$ Restando ambas expresiones, obtenemos que: $$f(x+h)-f(x-h) = 2f'(x)h+\color{red}{\dfrac{(f'''(\xi)+f'''(\eta))h^3}{6}}$$ $$\begin{array}{rll}\Rightarrow f'(x) & = & \dfrac{f(x+h)-f(x+h)}{2h} + \color{red}{\dfrac{(f'''(\xi)+f'''(\eta))h^2}{12}} \\ & = & \dfrac{f(x+h)-f(x+h)}{2h} + \color{red}{O(h^2)} \end{array}$$

Luego, $f'(x)$ puede aproximarse como: $$f'(x) \approx \dfrac{f(x+h)-f(x+h)}{2h}$$

Derivadas de funciones en varias variables

Derivadas parciales

Sean $f\colon\mathbb{R}^n\rightarrow \mathbb{R}$, $f\in C^1$, $h>0$, podemos usar las ideas proporcionadas por los métodos antes vistos para calcular las derivadas parciales. Veamos, por ejemplo, cómo aplicar el método de centradas para aproximar la derivada parcial con respecto a $x_i$: $$\dfrac{\partial f}{\partial x_i}(x) \approx \dfrac{f(x+he_i)-f(x-he_i)}{2h}$$

Gradiente

Una vez que tenemos las aproximaciones de las derivadas parciales, calcular el gradiente de $f$ es trivial, teniendo en cuenta que: $$\nabla f (x) = \left(\dfrac{\partial f}{\partial x_1}(x), \dots, \dfrac{\partial f}{\partial x_n}(x)\right) $$

Hessiano

El próximo paso, aproximar el valor del Hessiano en cierto $x$, también es relativamente sencillo, si tenemos en cuenta que: $$Hf(x)={\begin{pmatrix}{\frac {\partial ^{2}f}{\partial x_{1}^{2}}}&{\frac {\partial ^{2}f}{\partial x_{1}\partial x_{2}}}&\cdots &{\frac {\partial ^{2}f}{\partial x_{1}\partial x_{n}}}\\{\frac {\partial ^{2}f}{\partial x_{2}\partial x_{1}}}&{\frac {\partial ^{2}f}{\partial x_{2}^{2}}}&\cdots &{\frac {\partial ^{2}f}{\partial x_{2}\partial x_{n}}}\\\vdots &\vdots &\ddots &\vdots \\{\frac {\partial ^{2}f}{\partial x_{n}\partial x_{1}}}&{\frac {\partial ^{2}f}{\partial x_{n}\partial x_{2}}}&\cdots &{\frac {\partial ^{2}f}{\partial x_{n}^{2}}}\end{pmatrix}} = {\begin{pmatrix}{\nabla \frac{\partial f}{\partial x_1}}\\{\vdots}\\{{\nabla \frac{\partial f}{\partial x_n}}}\end{pmatrix}}$$

Observar que: $$\dfrac{\partial f}{\partial x_i \partial x_j}(x) = \dfrac{\partial}{\partial x_j}\left(\dfrac{\partial f}{\partial x_i}\right)(x) \approx \frac{\dfrac{\partial f}{\partial x_i}(x+he_j)-\dfrac{\partial f}{\partial x_i}(x-he_j)}{2h}$$

Algoritmo de Newton generalizado

Sean $f:\mathbb{R}^n\rightarrow\mathbb{R}$, $f\in C^3$, $x_0\in\mathbb{R}^n$ por Polinomio de Taylor centrado en $x_0$ se obtiene la siguiente aproximación:

$$f(x) \approx f(x_0)+\nabla f(x_0)(x-x_0)+\frac{1}{2}(x-x_0)^tHf(x_0)(x-x_0) = P_{x_0}(x)$$

Nos gustaría encontrar $x$ tal que minimice a $P_{x_0}$, entonces buscamos un $x$ que cumpla con las condiciones necesarias de primer orden:

$$0 = \nabla P_{x_0}(x) = \nabla f(x_0) + Hf(x_0)(x-x_0) \Leftrightarrow \color{red}{Hf(x_0)(x-x_0) = -\nabla f(x_0)} \Leftrightarrow x = x_0 - Hf(x_0)^{-1}\nabla f(x_0)$$

Observar que si $Hf(x_0)$ es definida positiva, podemos resolver el sistema en rojo utilizando descomposición de Cholesky o descomposición QR. Es decir, resolvemos $Hf(x_0)d = -\nabla f(x_0)$ y resulta entonces que $x = x_0+d$.

El algoritmo de Newton generalizado consiste en iterar la idea anterior, de ser posible, hasta alcanzar un mínimo: $$x_{n+1} = x_n - Hf(x_n)^{-1}\nabla f(x_n)$$

Criterios de parada

Como todo algoritmo iterativo, debemos establecer uno o varios criterios de parada. Por ejemplo:

  • $\lVert x_{n+1}-x_n \rVert < \epsilon $ (sucesión de Cauchy)

  • $\lVert \nabla f(x_n) \rVert < \epsilon $ (el gradiente está convergiendo a $0$)

  • $iter > M $ (Límite al número de iteraciones)

Python: funciones anónimas

En Python, el comando lambda nos permite definir funciones anónimas (es decir, funciones que no definimos utilizando def). Una de las ventajas que proporcionan es la posibilidad de definir nuevas funciones a partir de otras, pero manteniendo fijos ciertos argumentos. Por ejemplo, en nuestro código definimos la siguiente función:

In [1]:
def pot(base, potencia):
    return base**potencia

Supongamos ahora que en alguna parte del código necesitamos una nueva función que calcule sólo las potencias de $2$. Aquí puede resultar útil declarar una función anónima:

In [3]:
g = lambda x: pot(2,x)

Ahora g es una función. Notar que 2 está fijo, entonces g calcula $2^x$. Veámoslo:

In [4]:
print(g(2))
print(g(4))
print(g(-1))

4
16
0.5

Otra forma de pensarlo es la siguiente. Supongamos que definimos la funcion gradiente, que toma como input una función f y un vector x:

In [ ]:
def gradiente(func,x):
    # Implementacion de la aproximacion al gradiente

Si querríamos pensarlo conceptualmente, si llamamos $grad$ a la función que definimos, entonces

$$grad: Dif\times \mathbb{R}^n\rightarrow \mathbb{R}$$

donde $Dif$ es el conjunto funciones diferenciables. Ahora bien, si ahora definimos la siguiente función anónima:

In [ ]:
grad_f = lambda x: gradiente(f,x) # f ya debe estar definida en el codigo

efectivamente ahora grad_f$=\nabla f$

Una vez que uno implementa las funciones que calculan las derivadas parciales y el gradiente de una función, resulta útil recurrir a las funciones anónimas para implementar la aproximación del Hessiano: como la i-ésima fila del Hessiano es el gradiente de $\frac{\partial f}{\partial x_i}$, en la definición de la función hessiano podemos usar

In [ ]:
g = lambda x: derivada_parcial(f, x, i)

y utilizar la función gradiente sobre g para obtener la i-ésima fila del Hessiano.

Ejercicios

  • Implementar una función que aproxime la derivada de funciones en una variable, para cada uno de los métodos.
  • Implementar una función derivada_parcial que aproxime la i-ésima derivada parcial; debe recibir como input una función f, un vector a y la coordenada $i$ sobre la que se derivará parcialmente (¡recordar que la primera coordenada es $0$!). Debe devolver la aproximación de $\frac{\partial f}{\partial x_i}(a)$. Pueden aproximar con cualquiera de los métodos (o hacer una función para cada uno!)
  • Implementar la función gradiente. Debe recibir como input una función f y un vector a y devolver la aproximación de $\nabla f(x)$
  • Implementar la funcion hessiano. Debe recibir como imput una función f y un vector a y devolver la aproximación de $Hf(a)$
  • Implementar el algoritmo de Newton generalizado, que reciba una función f y un vector inicial x0 y devuelva el minimizador de la función.

Tips

  • Usar funciones anónimas!
  • (Opcional) En algoritmo de Newton generalizado: se puede implementar una función que devuelva si una matriz es definida positiva. Si lo es, resolver $Hf(x_n)d = -\nabla f(x_n)$ utilizando Cholesky y definir $x_{n+1} = x_n + d$.
  • (General) En algoritmo de Newton generalizado: resolver $Hf(x_n)d = -\nabla f(x_n)$ con el comando solve de numpy.linalg y definir $x_{n+1} = x_n + d$. Si uno quiere la solución de $Ax=b$, se utiliza solve de la siguiente manera x = solve(A,b).
In [ ]:
from numpy.linalg import solve
  • En el caso de la aproximación de las derivadas parciales, es interesante recalcularlas para distintos $h$ y en cierta manera asegurar una buena aproximación. Es decir, vamos achicando el $h$ hasta que los resultados sean muy parecidos (o iguales) o hasta que la división no esté definida
In [ ]:
from numpy.linalg import norm

def derivada_parcial(f,x,i):
    h = 0.1
    # bla ....
    z = #formula del metodo
    h = h/2
    y = #formula del metodo
    error = norm(y-z)
    eps = 1e-8
    while error>eps and (y != np.nan) and (y != np.inf) and y!= 0:
        error = norm(y-z)
        z = y
        h = h/2
        y = #formula formula del metodo
    return z