Revisión de Cálculo para Optimización

Mariano Rivera

%matplotlib inline

from sympy import *
init_printing()

x = Symbol('x')
y = Symbol('y')
z = Symbol('z')

Gradiente, Jacobiano y Hessiano

Derivada parcial

Calcular la derivada parcial de una función que depende de varias variables c.r.a (con respecto a) una variable, corresponde a calcular la derivada de la función respecto a esa variable asuminedo las demás constante. Esto es:
$\frac{\partial f(x_1, x_2, \ldots, x_i,\ldots, x_n)}{\partial x_i} \overset{def}{=} \lim_{h \rightarrow 0}\frac{ f(x_1, x_2, \ldots, x_i+h,\ldots, x_n) - f(x_1, x_2, \ldots, x_i,\ldots, x_n) }{h}.$

O mas compacto
$\boxed{ \frac{\partial f(x)}{\partial x_i} \overset{def}{=} \lim_{h \rightarrow 0}\frac{ f(x+h \,e_i)-f(x)}{h} }$
donde
$x = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ \vdots \\ x_n \\ \end{bmatrix}$

f = x**2 + 2*y**2 + 3*z**2 - 6*x**2*z + 4*y*x + 5*y*z**2 - 3*x*y*z
f

$- 6 x^{2} z + x^{2} - 3 x y z + 4 x y + 2 y^{2} + 5 y z^{2} + 3 z^{2}$

diff(f,x)

$- 12 x z + 2 x - 3 y z + 4 y$

diff(f,z)

$- 6 x^{2} - 3 x y + 10 y z + 6 z$

f = x**2 + 2*y**2 + 3*z**2
f

$x^{2} + 2 y^{2} + 3 z^{2}$

El operador Gradiente $\nabla$:
$\nabla f(x) \overset{def}{=} \begin{bmatrix} \frac{\partial f(x) }{ \partial x_1} \\ \frac{\partial f(x) }{ \partial x_2} \\ \frac{\partial f(x) }{ \partial x_3} \\ \vdots \\ \frac{\partial f(x) }{ \partial x_n} \\ \end{bmatrix}$

Podemos ver al gradiente, intuitivamente, que al operar sobre una función y la “expande” en un vector columna que contienen sus derivadas parciales.

Matriz Jacobiana o Jacobiano de la función $f$:

$J_f \overset{def}{=} \begin{bmatrix} \frac{\partial f(x) }{ \partial x_1}, & \frac{\partial f(x) }{ \partial x_2}, & \frac{\partial f(x) }{ \partial x_3}, & \ldots, & \frac{\partial f(x) }{ \partial x_n} \\ \end{bmatrix} \\$
$\;\;\;\;\;\;\;\;\;= \begin{bmatrix} \frac{\partial f_1(x) }{ \partial x_1} &, \frac{\partial f_1(x) }{ \partial x_2}, & \frac{\partial f_1(x) }{ \partial x_3}, &\ldots, &\frac{\partial f_1(x) }{ \partial x_n} \\ \frac{\partial f_2(x) }{ \partial x_1}, &\frac{\partial f_2(x) }{ \partial x_2}, & \frac{\partial f_2(x) }{ \partial x_3},& \ldots, &\frac{\partial f_2(x) }{ \partial x_n} \\ \frac{\partial f_3(x) }{ \partial x_1}, &\frac{\partial f_3(x) }{ \partial x_2}, & \frac{\partial f_3(x) }{ \partial x_3}, &\ldots, & \frac{\partial f_3(x) }{ \partial x_n} \\ \vdots, & \vdots& \vdots & \ddots & \vdots \\ \frac{\partial f_m(x) }{ \partial x_1}, &\frac{\partial f_m(x) }{ \partial x_2}, & \frac{\partial f_m(x) }{ \partial x_3}, &\ldots, & \frac{\partial f_m(x) }{ \partial x_n} \end{bmatrix}$

Podemos ver intuitivamente el Jacobiano como un operador que “expande” cada renglón de la función vectorial $f$ en sus derivadas parciales.

Luego si $f \equiv f_1$:
$\nabla f = J_f^\top$

sbl = Matrix([x,y,z])
sbl

$\left[\begin{matrix}x\\y\\z\end{matrix}\right]$

f = Matrix([f])
f

$\left[\begin{matrix}x^{2} + 2 y^{2} + 3 z^{2}\end{matrix}\right]$

G=f.jacobian(sbl).T
G

$\left[\begin{matrix}2 x\\4 y\\6 z\end{matrix}\right]$

Recordando la definición de Gradiente de $f$
$\nabla f(x) = \begin{bmatrix} \frac{\partial f(x) }{ \partial x_1} \\ \frac{\partial f(x) }{ \partial x_2} \\ \frac{\partial f(x) }{ \partial x_3} \\ \vdots \\ \frac{\partial f(x) }{ \partial x_n} \\ \end{bmatrix}$

Matrix Hessiana o Hessiano de $f$:

• El Hessiano es una matriz de segundas derivadas

• Es simétrico si
  $\frac{\partial^2 }{\partial x_i x_j} f = \frac{\partial^2 }{\partial x_j x_i} f,$
  esto es $f$ es 2 veces continuamente diferenciable $f \in \mathcal{C}^2$.

Note que el Hessiano, se puede ver como aplicar el operador Jacobiano al Gradiente:
$\boxed{ H_f \overset{def}{=} J \{ \nabla f\} = J\{ J_f^\top\} }$

Sobre la Notación:

$J\{ \cdot \}$ se refiere al operador

$J_f$ se refiere a la Matrix

H = G.jacobian(sbl)
H

$\left[\begin{matrix}2 & 0 & 0\\0 & 4 & 0\\0 & 0 & 6\end{matrix}\right]$

det(H)

$48$

Probablemente, nos sintamos mas cómodos calculando el gradiente con un operador especifico de Gradiente que con el de Jacobiano. Para tal caso, definimos lo siguiente

def Grad(f, variables):
    return Matrix([f]).jacobian(Matrix(symbols(variables))).T

La divergencia se define como
$\nabla \cdot f \overset{def}{=} \mathbf{1}^\top \nabla f$

def Div(f, variables):
     return sum(Grad(f,variables))

Hessiano se define como:
$\nabla^2 f \overset{def}{=} \nabla^\top \nabla f$
Como en python no podemos indicarle al operador matrix.jacobian como haga la expansion. Nuestro Grad, siempre calculará derivadas por columnas asi que jugamos un poco con el orden de los operadores para implementrar el Hessiano.

Si $f$ es dos veces continua y diferenciable, entonces tenemos
$\nabla^2 f = ( \nabla^2 f)^\top = (\nabla^\top \nabla f)^\top = \nabla \nabla f^\top$

def Hessian(f, variables):
    g = Grad(f,variables)
    return Grad(g.T, variables)

Probemos nuestras implementaciones

from sympy import *
init_printing()

x1 = Symbol('x1')
x2 = Symbol('x2')

f = (x1 - x2)**2+ x1**3
f

$x_{1}^{3} + \left(x_{1} - x_{2}\right)^{2}$

Grad(f,'x1, x2')

$\left[\begin{matrix}3 x_{1}^{2} + 2 x_{1} - 2 x_{2}\\- 2 x_{1} + 2 x_{2}\end{matrix}\right]$

Div(f,'x1 x2')

$3 x_{1}^{2}$

Hessian(f, 'x1 x2')

$\left[\begin{matrix}6 x_{1} + 2 & -2\\-2 & 2\end{matrix}\right]$

Revisando propiedades de las derivadas que involucran Algebra Lineal

$\nabla_x y^\top x$

from sympy import *
init_printing()
from sympy import Matrix, latex, symbols

A = Matrix(symbols('a_1:4\,1:4')).reshape(3,3)
x = Matrix(symbols('x_1:4')).reshape(3,1)
y = Matrix(symbols('y_1:4')).reshape(3,1)

(y.T*x).jacobian(x).T

$\left[\begin{matrix}y_{1}\\y_{2}\\y_{3}\end{matrix}\right]$

(A*x).jacobian(x).T

$\left[\begin{matrix}a_{1,1} & a_{2,1} & a_{3,1}\\a_{1,2} & a_{2,2} & a_{3,2}\\a_{1,3} & a_{2,3} & a_{3,3}\end{matrix}\right]$

(x.T*A*x).jacobian(x).T

$\left[\begin{matrix}2 a_{1,1} x_{1} + a_{1,2} x_{2} + a_{1,3} x_{3} + a_{2,1} x_{2} + a_{3,1} x_{3}\\a_{1,2} x_{1} + a_{2,1} x_{1} + 2 a_{2,2} x_{2} + a_{2,3} x_{3} + a_{3,2} x_{3}\\a_{1,3} x_{1} + a_{2,3} x_{2} + a_{3,1} x_{1} + a_{3,2} x_{2} + 2 a_{3,3} x_{3}\end{matrix}\right]$

A*x+A.T*x

$\left[\begin{matrix}2 a_{1,1} x_{1} + a_{1,2} x_{2} + a_{1,3} x_{3} + a_{2,1} x_{2} + a_{3,1} x_{3}\\a_{1,2} x_{1} + a_{2,1} x_{1} + 2 a_{2,2} x_{2} + a_{2,3} x_{3} + a_{3,2} x_{3}\\a_{1,3} x_{1} + a_{2,3} x_{2} + a_{3,1} x_{1} + a_{3,2} x_{2} + 2 a_{3,3} x_{3}\end{matrix}\right]$

(x.T*A*x).jacobian(x).T- A*x-A.T*x

$\left[\begin{matrix}0\\0\\0\end{matrix}\right]$

J= (x.T*A*x).jacobian(x)
J.jacobian(x)

$\left[\begin{matrix}2 a_{1,1} & a_{1,2} + a_{2,1} & a_{1,3} + a_{3,1}\\a_{1,2} + a_{2,1} & 2 a_{2,2} & a_{2,3} + a_{3,2}\\a_{1,3} + a_{3,1} & a_{2,3} + a_{3,2} & 2 a_{3,3}\end{matrix}\right]$

A+A.T

$\left[\begin{matrix}2 a_{1,1} & a_{1,2} + a_{2,1} & a_{1,3} + a_{3,1}\\a_{1,2} + a_{2,1} & 2 a_{2,2} & a_{2,3} + a_{3,2}\\a_{1,3} + a_{3,1} & a_{2,3} + a_{3,2} & 2 a_{3,3}\end{matrix}\right]$

Optimización Convexa

Conjunto Convexo

Definiciones

Conjunto Convexo. $\Omega \subseteq \mathbb{R}^n$ es un conjunto convexo si
$\alpha x + (1-\alpha) y \in \Omega, \;\;\; \forall x, y \in \Omega \;\; y\;\; \forall \alpha \in [0,1]$

Función Convexa. $f:\Omega \rightarrow \mathbb{R}$ es una función convexa si
$\alpha f(x) + (1-\alpha) f(y) \ge f ( \alpha x + (1-\alpha) y ), , \;\;\; \forall x, y \in \Omega$
y el dominio $\Omega$ es un conjunto convexo

$f(x) = |x|$

Función Estrictamente Convexa. $f:\Omega \rightarrow \mathbb{R}$ es una función convexa si
$\alpha f(x) + (1-\alpha) f(y) > f ( \alpha x + (1-\alpha) y ), , \;\;\; \forall x, y \in \Omega$
y el dominio $\Omega$ es un conjunto convexo

$f(x) = x^2$

Función que es una aproximación diferenciable al valor absoluto. Esta función aproxima al valor absoluto para $x >> \beta$
$f(x) = \sqrt{x^2+\beta^2}$

Optimalidad

Sea el problema de optimización
$\boxed{ \min_{x} \; f(x) }$

con $x\in \mathbb{R}^n$ y $f \in \mathcal{C}^2$.

Óptimo Global Estricto. $x^\ast \in \mathbb{R}^n$ es el óptimo global estricto si $f(x^\ast) < f(y), \;\; \forall y \in \mathbb{R}^n$.

Óptimo Global. $x^\ast \in \mathbb{R}^n$ es el óptimo global si $f(x^\ast) \le f(y), \;\; \forall y \in \mathbb{R}^n$.

Óptimo Local Estricto. $x^\ast \in \mathbb{R}^n$ es el óptimo local estricto si $f(x^\ast) < f(y), \;\; \forall y \in \mathcal{N}_{x^\ast}$; donde $\mathcal{N}_{x^\ast}$ es una vecindad abierta que contine a $x^\ast$.

Óptimo Local o mínimo local (m.l). $x^\ast \in \mathbb{R}^n$ es el óptimo local si $f(x^\ast) \le f(y), \;\; \forall y \in \mathcal{N}_{x^\ast}$; donde $\mathcal{N}_{x^\ast}$ es una vecindad abierta que contine a $x^\ast$.

Cuando implementamos algoritmos de optimización generalmente nos conformamos con encontrar un óptimo local como solución .

Dirección de Descenso

Dirección de descenso. $d$ es una dirección de descenso en el punto $x'$ para el problema propuestonde optimización si partiendo de x’ en la dirección d podemos mejorar el valor de la función objetivo, esto es:
$f(x'+\epsilon d) < f(x')$
para una $\epsilon$ suficientemente pequeña (sfmp).

Para encontrar que condición cumple $d$, usamos expansión en serie de Taylor:
$f(x') > f(x'+\epsilon d) = f(x') + \epsilon \nabla f(x') ^\top d + \epsilon^2 d^\top \nabla^2 f \, d + \ldots \\ \;\;\;\;\;\; \approx f(x') + \epsilon \nabla f(x') ^\top d.$
Dado que $\epsilon$ es pequeño: $| \epsilon |< | \epsilon^n |$ para $n>1$.
Luego restamos a ambos lados $f(x')$ y tenemos que:

La dirección de descenso $d$ debe cumplir con
$\boxed{ \nabla f(x')^ \top d < 0 }$

La línea punteada es tangente a la curva de nivel correspondiente a $f$ en $x'$

El gradiente $\nabla f(x')$ es normal a la curva de nivel y apunta en la dirección en que la función $f$ crece.

Gráficamente, $d$ debe estar en semiespacio opuesto a $\nabla f$ que es delimitado por la tangente a la curva de nivel de $f$ en $x'$.

Condición de Optimalidad de Primer Orden

Asumimos el problema de optimización sin restricciones de una función suave.

Una solución (m.l.) $x^\ast$ al problema de optimización es un punto en el cual no existe una dirección de descenso. Esto es:
$\nexists d : \nabla f(x^\ast)^ \top d < 0$
Esto solo ocurre si:
$\nabla f(x^\ast) = 0$

Condicion de optimalidad de primer orden. Si $x^\ast$ es un m.l. entonces $\nabla f(x^\ast) = 0$.

Si la función $f$ es convexa, entonces el m.l. es mínimo global y el punto $x^*$ que satisfaga $\nabla f(x^\ast) = 0$ es solución global al problema de optimización

En los máximos también el gradiente se hace cero

Pensemos que queremos encontrar una dirección de ascenso $b$ a partir del punto actual $x'$. Entonces $b$ debe cumplit con:
$f(x'+\epsilon b) > f(x')$
para una $\epsilon$ suficientemente pequeña (sfmp).

Usamos expanción en serie de Taylor de 1er orden:
$f(x') < f(x'+\epsilon d) \approx f(x') + \epsilon \nabla f(x') ^\top d.$
Como antes, restamos a ambos lados $f(x')$ y tenemos que una dirección de ascenso $b$ debe cumplir con
$\boxed{ \nabla f(x')^ \top b > 0 }$

Si x’ es un maximo local, entonces si existe una dirección de ascenso posible. La única forma de no entrar en contradicción con condición anterior, es que $\nabla f(x')=0$.

Condiciones de Optimalidad de Segundo Orden

En general (funciones no-convexas), la condición de optimalidad 1er orden es necesaria, pero no es suficiente para encontrar una solución.

Es decir, en un m.l. el gradiente se desvanece; pero también el gradiente se desvanece en otros puntos (máximos) que no son m.l.'s.

Para distinguir estos casos necesitamos la condiciones de segundo orden.

Condiciones de segundo orden sin restricciones

Sea el problema de optimización SIN restricciones

$\min_x \;\; f(x)$

La solución $x^*$ satisface la condición de optimalidad primer orden
$\nabla_x f(x^*) = {\bf 0}$

Esto evita que un método de descenso actualice el punto óptimo:
$x^\ast \leftarrow x^\ast - \alpha \nabla_x f(x^*)$

Sin embargo, usiando la approximación de series de Taylor de segundo orden notamos que para $\epsilon$ suficientemente pequeña (sfmp):
$f(x^*+ \epsilon p) = f(x^*) + \epsilon p^\top \nabla f(x^*) +\frac{1}{2} \epsilon^2\,p^\top \nabla^2 f(x^*) \, p \\ = f(x^*) +\frac{1}{2} \epsilon^2\,p^\top \nabla^2 f(x^*) \, p$
luego, dado que $x^*$ es óptimo:
$p^\top \nabla^2 f(x^*) \, p = \frac{2}{\epsilon^2} \left[ f(x^*+ \epsilon p) - f(x^*) \right] \ge 0, \;\; \forall p.$

Condición necesaria de segundo orden

TMA.
Si $x^*$ es solución del problema
$\min_x \;\; f(x)$
Entonces, $\nabla f(x^*)=0$ y el Hessiano $\nabla^2 f(x^*)$ en dicho óptimo es semidefinido positivo.

Proof.
De lo anteriormente visto tenemos $p^\top \nabla^2 f(x^*) \, p \ge 0, \;\; \forall p$.

Condición suficiente de segundo orden

TMA. Sea el problema
$\min_x \;\; f(x)$
con f suave (dos veces continuamente diferenciable)y $x'$ un punto en el cual

• el gradiente $\nabla f(x')=0$ y
• el Hessiano $\nabla^2 f(x')$ es definido positivo;

entonces, $x'$ es un m.l. del problema de optimización.

Proof. Sea el punto $x'$ tal que $\nabla f(x') =0$ y $\, p^\top \nabla^2 f(x') p >0$. Luego elegimos una bola abierta de radio $r$ al rededor de $x'$
$D = \{z \,:\, \|x'-z \|< r \},$
tal que $\nabla^2 f(x')$ se mantenga definido positivo. Tomado un vector cualquier vector $p$ tal que
$x'+ \epsilon p \in D$

Luego, para $\epsilon$ sfmp:
$f(x'+ \epsilon\,p) = f(x') + \epsilon\, p^\top \nabla f(x') + \frac{\epsilon^2}{2} p^\top \nabla^2 f(x') p$
Entonces, dado que:

• $\nabla f(x')=0\,$, y
• $\epsilon^2\, p^\top \nabla^2 f(x') p >0$;

tenemos que
$f(x'+\epsilon p) > f(x'),$
por lo que $x'$ es un m.l.

import numpy as np
import matplotlib.pyplot as plt

N=10
x= (np.arange(0,N)/N)**2
z= np.zeros(int(N/2))
f1 = np.concatenate([ x[::-1],z,x])
f2 = np.concatenate([-x[::-1],z,x])
f3 = np.concatenate([x[::-1],x])

fig =plt.figure(figsize=(15,5))

ax = fig.add_subplot(1,3,1)
ax.plot(f3,'b' )
ax.annotate(r'$\nabla f=0, \;\; \nabla^2f \succ 0$',
            xy=(9.5,0), xytext=(6.4,.2),
            arrowprops=dict(facecolor='red', shrink=1),
            )
ax = fig.add_subplot(1,3,2)
ax.plot(f1,'b' )
ax.annotate(r'$\nabla f=0, \;\; \nabla^2f = 0$', 
            xy=(12,0), xytext=(8,.15),
            arrowprops=dict(facecolor='red', shrink=1),
            )
ax = fig.add_subplot(1,3,3)
ax.plot(f2,'b' )
ax.annotate(r'$\nabla f=0, \;\; \nabla^2f = 0$',
            xy=(12,0), xytext=(8,.3),
            arrowprops=dict(facecolor='red', shrink=1),
            )

De los tres casos mostrados en la figura anterior, los dos promeros (de la izquierda) corresponden a mínimos locales (soluciones) y el tercero a un punto de mesa (no-solución).

La tercera gráfica (a la derecha) ilustra el porqué para evitar esta caso, en la condición suficiente de segundo orden tuvimos que asumir de inicio que el Hessiano era definido positivo en $x'$.

Lo que concluimos que:

• si $x^*$ es óptimo, entonces cumple $\nabla f^*=0, \;\; \nabla^2 f^* \succeq 0$ (condiciones necesarias)

• si se cumple $\nabla f^*=0, \;\; \nabla^2 f^* \succ 0$, entonces $x^*$ es óptimo (condiciones suficientes)