Notas para contenedor de docker:
Comando de docker para ejecución de la nota de forma local:
nota: cambiar <ruta a mi directorio> por la ruta de directorio que se desea mapear a /datos dentro del contenedor de docker.
docker run --rm -v <ruta a mi directorio>:/datos --name jupyterlab_numerical -p 8888:8888 -d palmoreck/jupyterlab_numerical:1.1.0password para jupyterlab: qwerty
Detener el contenedor de docker:
docker stop jupyterlab_numericalDocumentación de la imagen de docker palmoreck/jupyterlab_numerical:1.1.0 en liga.
Optimización de código ¿es optimización numérica?
Hasta este módulo hemos invertido buena parte del tiempo del curso en la eficiente implementación en el hardware que poseemos. Revisamos lo que estudia el análisis numérico o cómputo científico, definiciones de sistema de punto flotante, funciones, derivadas, integrales y métodos o algoritmos numéricos para su aproximación. Consideramos bottlenecks que pueden surgir en la implementación de los métodos o algoritmos y revisamos posibles opciones para encontrarlos y minimizarlos (optimización de código). Lo anterior lo resumimos con el uso de herramientas como: perfilamiento, integración de R y Python con C++ o C, cómputo en paralelo y uso del caché de forma eficiente al usar niveles altos en operaciones de BLAS (módulo I: cómputo científico y análisis numérico, módulo II: cómputo en paralelo, módulo III: cómputo matricial).
La optimización numérica no es optimización de código sin embargo se apoya enormemente de ella para la implementación de sus métodos o algoritmos en la(s) máquina(s) para resolver problemas que surgen en tal rama de las matemáticas aplicadas. A la implementación y simulación en el desarrollo de los métodos o algoritmos del análisis numérico o cómputo científico típicamente se le acompaña de estudios que realizan benchmarks) y perfilamiento (mediciones de tiempo y memoria, por ejemplo) con el objetivo de tener software confiable y eficiente en la práctica. Esto lo encontramos también en la rama de optimización numérica con los métodos o algoritmos que son desarrollados e implementados.
Métodos o algoritmos numéricos en big data
La implementación de los métodos o algoritmos en el contexto de grandes cantidades de datos o big data es crítica al ir a la práctica pues de esto depende que nuestra(s) máquina(s) tarde meses, semanas, días u horas para resolver problemas que se presentan en este contexto. La ciencia de datos apunta al desarrollo de técnicas y se apoya de aplicaciones de machine learning para la extracción de conocimiento útil y toma como fuente de información las grandes cantidades de datos.
Una gran cantidad de aplicaciones plantean problemas de optimización. Tenemos problemas básicos que se presentan en cursos iniciales de cálculo:
Una caja con base y tapa cuadradas debe tener un volumen de $100 cm^3$. Encuentre las dimensiones de la caja que minimicen la cantidad de material.
Y tenemos más especializados que encontramos en áreas como estadística, ingeniería, finanzas o machine learning:
Ajustar un modelo de regresión lineal a un conjunto de datos.
Buscar la mejor forma de invertir un capital en un conjunto de activos.
Elección del ancho y largo de un dispositivo en un circuito electrónico.
Ajustar un modelo que clasifique un conjunto de datos.
En general un problema de optimización matemática o numérica tiene la forma:
$$\displaystyle \min f_o(x)$$donde: $x=(x_1,x_2,\dots, x_n)^T$ es la variable de optimización del problema, la función $f_o: \mathbb{R}^{n} \rightarrow \mathbb{R}$ es la función objetivo, las funciones $f_i: \mathbb{R}^n \rightarrow \mathbb{R}, i=1,\dots,m$ son las funciones de restricción (aquí se colocan únicamente desigualdades pero pueden ser sólo igualdades o bien una combinación de ellas) y las constantes $b_1,b_2,\dots, b_m$ son los límites o cotas de las restricciones.
Un vector $x^* \in \mathbb{R}^n$ es nombrado óptimo o solución del problema anterior si tiene el valor más pequeño de entre todos los vectores $x \in \mathbb{R}^n$ que satisfacen las restricciones. Por ejemplo, si $z \in \mathbb{R}^n$ satisface $f_1(z) \leq b_1, f_2(z) \leq b_2, \dots, f_m(z) \leq b_m$ y $x^*$ es óptimo entonces $f_o(z) \geq f_o(x^*)$.
Comentarios:
En el módulo IV del curso revisaremos métodos o algoritmos de optimización para funciones objetivo $f_o: \mathbb{R}^n \rightarrow \mathbb{R}$. Sin embargo, hay formulaciones que utilizan $f_o: \mathbb{R}^n \rightarrow \mathbb{R}^q$. Tales formulaciones pueden hallarlas en la optimización multicriterio, multiobjetivo, vectorial o también nombrada Pareto, ver Multi objective optimization.
Obsérvese que el problema de optimización definido utiliza una forma de minimización y no de maximización. Típicamente en la literatura por convención se consideran problemas de este tipo. Además minimizar $f_o$ y maximizar $-f_o$ son problemas de optimización equivalentes*.
*A grandes rasgos dos problemas de optimización son equivalentes si con la solución de uno de ellos se obtiene la solución del otro y viceversa.
Ejemplo:
1) $$\displaystyle \min_{x \in \mathbb{R}^n} ||x||_2$$
con $A \in \mathbb{R}^{m \times n}, b \in \mathbb{R}^m$. En este problema buscamos el vector $x$ que es solución del problema $Ax \leq b$ con mínima norma Euclidiana. La función objetivo es $f_o(x)=||x||_2$, las funciones de restricción son las desigualdades lineales $f_i(x) = a_i^Tx \leq b_i$ con $a_i$ $i$-ésimo renglón de $A$ y $b_i$ $i$-ésima componente de $b$, $\forall i=1,\dots,m$.
Comentario: un problema similar al anterior lo podemos encontrar en resolver el sistema de ecuaciones lineales $Ax=b$ underdetermined en el que $m < n$ y se busca el vector $x$ con mínima norma Euclidiana que satisfaga tal sistema. Tal sistema puede tener infinitas soluciones o ninguna solución, ver 3.3.Solucion_de_SEL_y_FM.
2) Encuentra el punto en la gráfica de $y=x^2$ que es más cercano al punto $P=(1,0)$ bajo la norma Euclidiana.
Deseamos minimizar la cantidad $||(1,0)-(x,y)||_2$. Además $y = y(x)$ por lo que reescribiendo lo anterior se tiene $f_o(x) = ||(1,0)-(x,x^2)||_2=||(1-x,-x^2)||_2=\sqrt{(1-x)^2+x^4}$. Entonces el problema de optimización (sin restricciones) es:
$$\displaystyle \min_{x \in \text{dom}f_o}\sqrt{(1-x)^2+x^4}$$En esta sección relacionamos a machine learning con la optimización y se describen diferentes enfoques que se han propuesto para aplicaciones de machine learning con métodos de optimización. Lo siguiente no pretende ser una exposición extensa ni completa sobre machine learning, ustedes llevan materias que se enfocan esencialmente a definir esta área, sus objetivos y conceptos más importantes.
En la ciencia de datos se utilizan las aplicaciones desarrolladas en machine learning por ejemplo:
Clasificación de documentos o textos: detección de spam.
Procesamiento de lenguaje natural: named-entity recognition.
Visión por computadora: reconocimiento de rostros o imágenes.
Detección de fraude.
Diagnóstico médico.
Las aplicaciones anteriores involucran problemas como son:
Clasificación.
Regresión.
Ranking.
Clustering.
Reducción de la dimensionalidad.
En cada una de las aplicaciones o problemas anteriores se utilizan funciones de pérdida que guían el proceso de aprendizaje. Tal proceso involucra optimización parámetros de la función de pérdida. Por ejemplo, si la función de pérdida en un problema de regresión es una pérdida cuadrática $\mathcal{L}(y,\hat{y}) = (\hat{y}-y)^2$ con $\hat{y} = \hat{\beta}_0 + \hat{\beta_1}x$, entonces el vector de parámetros a optimizar (aprender) es $ \beta= \left[ \begin{array}{c} \beta_0\\ \beta_1 \end{array} \right] $.
Machine learning no sólo se apoya de la optimización pues es un área de Inteligencia Artificial* que utiliza técnicas estadísticas para el diseño de sistemas capaces de aplicaciones como las escritas anteriormente, de modo que hoy en día tenemos statistical machine learning. No obstante, uno de los pilares de machine learning o statistical machine learning es la optimización.
*La IA o inteligencia artificial es una rama de las ciencias de la computación que atrajo un gran interés en $1950$.
Machine learning o statistical machine learning se apoya de las formulaciones y algoritmos en optimización. Sin embargo, también ha contribuido a ésta área desarrollando nuevos enfoques en los métodos o algoritmos para el tratamiento de grandes cantidades de datos o big data y estableciendo retos significativos no presentes en problemas clásicos de optimización. De hecho, al revisar literatura que intersecta estas dos disciplinas encontramos comunidades científicas que desarrollan o utilizan métodos o algoritmos exactos (ver Exact algorithm) y otras que utilizan métodos de optimización estocástica (ver Stochastic optimization y Stochastic approximation) basados en métodos o algoritmos aproximados (ver Approximation algorithm). Hoy en día es común encontrar estudios que hacen referencia a modelos o métodos de aprendizaje.
Como ejemplo de lo anterior considérese la técnica de regularización que en machine learning se utiliza para encontrar soluciones que generalicen y provean una explicación no compleja del fenómeno en estudio. La regularización sigue el principio de la navaja de Occam, ver Occam's razor: para cualquier conjunto de observaciones en general se prefieren explicaciones simples a explicaciones más complicadas. Aunque la técnica de regularización es conocida en optimización, han sido varias las aplicaciones de machine learning las que la han posicionado como clave.
El inicio del siglo XXI estuvo marcado, entre otros temas, por un incremento significativo en la generación de información. Esto puede contrastarse con el desarrollo de los procesadores de las máquinas, que como se revisó en el módulo II del curso en el tema 2.1.Un_poco_de_historia_y_generalidades, tuvo un menor performance al del siglo XX. Asimismo, las mejoras en dispositivos de almacenamiento o storage y sistemas de networking abarató costos de almacenamiento y permitió tal incremento de información. En este contexto, los modelos y métodos de statistical machine learning se vieron limitados por el tiempo de cómputo y no por el tamaño de muestra. La conclusión de esto fue una inclinación en la comunidad científica por el diseño o uso de métodos o modelos para procesar grandes cantidades de datos usando recursos computacionales comparativamente menores.
Un ejemplo de lo anterior se observa en métodos de optimización desarrollados en la década de los $50$'s. Mientras que métodos tradicionales en optimización basados en el cálculo del gradiente y la Hessiana de una función son efectivos para problemas de aprendizaje small-scale (en los que utilizamos un enfoque en batch o por lote), en el contexto del aprendizaje large-scale, el método de gradiente estocástico* se posicionó en el centro de discusiones a inicios del siglo XXI.
* El método de gradiente estocástico fue propuesto por Robbins y Monro en 1951, es un algoritmo estocástico. Ver Stochastic gradient descent.
Tradicionalmente en optimización, la búsqueda del (o los) óptimo(s) involucran el cálculo de información de primer o segundo orden (ver 1.4.Polinomios_de_Taylor_y_diferenciacion_numerica) de la función $f_o$:
Ejemplo:
1) Calcular $\nabla f(x), \nabla^2f(x)$ con $f: \mathbb{R}^4 \rightarrow \mathbb{R}$, dada por $f(x) = (x_1-2)^2+(2-x_2)^2+x_3^2+x_4^4$ en el punto $x_0=(1.5,1.5,1.5,1.5)^T$.
Solución:
La información de primer y segundo orden la constituyen el gradiente de $f$, $\nabla f(x)$, y la matriz Hessiana de $f$, $\nabla^2f(x)$. Obsérvese que el almacenamiento de la Hessiana involucra $\mathcal{O}(n^2)$ entradas. En los métodos clásicos de optimización se utilizan el gradiente y la Hessiana para encontrar el mínimo de funciones. La Hessiana se utiliza para resolver un sistema de ecuaciones lineales asociado.
2) Encontrar el mínimo de $f$.
Solución:
In [1]:
import numpy as np
In [2]:
x0=np.array([1.5,1.5,1.5,1.5])
In [3]:
gf= lambda x: np.array([2*(x[0]-2),
-2*(2-x[1]),
2*x[2],
4*x[3]**3])
In [4]:
Hf = lambda x: np.array([[2, 0, 0 ,0],
[0, 2, 0, 0],
[0, 0, 2, 0],
[0, 0, 0, 12*x[2]**2]])
In [5]:
gf(x0)
Out[5]:
In [6]:
Hf(x0)
Out[6]:
Como $f$ es una función convexa (definida más adelante) se tiene que su óptimo se obtiene igualando y resolviendo la ecuación no lineal $\nabla f(x) = 0$ :
$$\nabla f(x) = \left[ \begin{array}{c} 2(x_1-2) \\ -2(2-x_2)\\ 2x_3\\ 4x_4^3 \end{array} \right] = 0 $$El óptimo $x^* \in \mathbb{R}^4$ está dado por:
$$x^*= \left[ \begin{array}{c} 2\\ 2\\ 0\\ 0 \end{array} \right] $$¿Cómo encontramos numéricamente el óptimo?
Forma1: con el gradiente de $f$
Numéricamente se puede utilizar un método iterativo en el que iniciamos con un punto inicial $x^{(0)}$ y las actualizaciones las realizamos con el gradiente:
$$x^{(k+1)} = x^{(k)} - \nabla f(x^{(k)})$$para $k=0,1,2,\dots,$.
Obs: A la iteración anterior se le añade la búsqueda de línea por backtracking, ver Backtracking line search para determinar el máximo paso a realizar en la dirección de descenso (que en este caso es el gradiente).
En el ejemplo tomando $x^{(0)} = (5,5,1,0)^T$ se tiene:
In [7]:
x_0 = np.array([5,5,1,0])
In [8]:
x_1 = x_0 - gf(x_0)
In [9]:
x_1
Out[9]:
In [10]:
x_2 = x_1 - gf(x_1)
In [11]:
x_2
Out[11]:
In [12]:
x_3 = x_2 - gf(x_2)
In [13]:
x_3
Out[13]:
In [14]:
x_4 = x_3 - gf(x_3)
In [15]:
x_4
Out[15]:
y aquí nos quedaremos ciclando hasta el infinito...
Forma2: con el gradiente y la Hessiana de $f$
Otra opción es utilizar la información de segundo orden con la Hessiana y considerar una actualización:
para $k=0,1,2,\dots,$.
Obs: A la iteración anterior se le añade la búsqueda de línea por backtracking, ver Backtracking line search para determinar el máximo paso a realizar en la dirección de descenso (que en este caso es la dirección de Newton).
En el ejemplo tomando $x^{(0)} = (5,5,1,0)^T$ se tiene lo siguiente. Recuérdese que no invertimos a la matriz y en su lugar resolvemos un sistema de ecuaciones lineales:
In [16]:
x_0 = np.array([5,5,1,0])
In [17]:
x_1 = x_0 - np.linalg.solve(Hf(x_0),gf(x_0))
In [18]:
x_1
Out[18]:
Comentarios: de acuerdo al ejemplo anterior:
Utilizar información de primer o segundo orden nos ayuda a encontrar óptimo(s) de funciones.
Encontrar al óptimo involucró un método iterativo.
En términos coloquiales y de forma simplificada, una dirección de descenso es aquella que al moverse de un punto a otro en tal dirección, el valor de $f_o$ decrece:
En el dibujo anterior $\hat{f}$ es un modelo cuadrático, $\Delta x_{nt}$ es dirección de descenso de Newton y $x^*$ es el óptimo de $f$. Del punto $(x,f(x))$ nos debemos mover al punto $(x+\Delta x_{nt}, f(x + \Delta x_{nt}))$ para llegar al óptimo y el valor de $f$ decrece: $f(x+\Delta x_{nt}) < f(x)$.
Con la información de primer orden no alcanzamos al óptimo (de hecho se cicla el método iterativo propuesto) pero con la de segundo orden sí lo alcanzamos en una iteración y tuvimos que resolver un sistema de ecuaciones lineales.
Si consideramos una reescritura en la actualización de descenso en gradiente:
de la forma:
para $t_{k} > 0 $. Con $t_0=0.5$ llegamos al óptimo en una iteración:
In [19]:
t_0=0.5
Comentario: más adelante veremos cómo obtener las $t_{k}$s de la búsqueda de línea por backtracking
In [20]:
gf(x_0)
Out[20]:
In [21]:
x_1 = x_0 - t_0*gf(x_0)
In [22]:
x_1
Out[22]:
para $k=0,1,2,\dots$ y $t_{k} >0$.
Ejemplo: regresión lineal
comentario: para este ejemplo la variable de optimización no será $x$, será $\beta$.
Supóngase que se han realizado mediciones de un fenómeno de interés en diferentes puntos $x_i$'s resultando en cantidades $y_i$'s $\forall i=0,1,\dots, m$ (se tienen $m+1$ puntos) y además las $y_i$'s contienen un ruido aleatorio causado por errores de medición:
El objetivo de los mínimos cuadrados es construir una curva, $f(x|\beta)$ que "mejor" se ajuste a los datos $(x_i,y_i)$, $\forall i=0,1,\dots,m$. El término de "mejor" se refiere a que la suma: $$\displaystyle \sum_{i=0}^m (y_i -f(x_i|\beta))^2$$ sea lo más pequeña posible, esto es, a que la suma de las distancias verticales entre $y_i$ y $f(x_i|\beta)$ $\forall i=0,1,\dots,m$ al cuadrado sea mínima. Por ejemplo:
Comentarios:
La notación $f(x|\beta)$ se utiliza para denotar que $\beta$ es un vector de parámetros a estimar, en específico $\beta_0, \beta_1, \dots \beta_n$, esto es: $n+1$ parámetros a estimar.
La variable de optimización es $\beta$.
donde: $r(\beta) = y-A\beta$, el vector $y \in \mathbb{R}^m$ contiene las entradas $y_i$'s y la matriz $A \in \mathbb{R}^{m \times n}$ contiene a las entradas $x_i$'s o funciones de éstas $\forall i=0,1,\dots,m$. Por el dibujo se tiene que cumplir que $A^Tr(\beta)=0$, esto es: las columnas de $A$ son ortogonales a $r(\beta)$. La condición anterior conduce a las ecuaciones normales:
$$0=A^Tr(\beta)=A^T(y-A\beta)=A^Ty-A^TA\beta.$$En los mínimos cuadrados lineales se supone: $f(x|\beta) = \displaystyle \sum_{j=0}^n\beta_j\phi_j(x)$ con $\phi_j: \mathbb{R} \rightarrow \mathbb{R}$ funciones conocidas por lo que se tiene una gran flexibilidad para el proceso de ajuste. Con las funciones $\phi_j (\cdot)$ se construye a la matriz $A$.
Obs:
En lo siguiente se asume $n+1 \leq m+1$ (tenemos más puntos $(x_i,y_i)$'s que parámetros a estimar).
Para realizar el ajuste de mínimos cuadrados lineales se utilizan las ecuaciones normales: $$A^TA\beta=A^Ty$$ donde: $A$ se construye con las $\phi_j$'s evaluadas en los puntos $x_i$'s, el vector $\beta$ contiene a los parámetros $\beta_j$'s a estimar y el vector $y$, la variable respuesta, se construye con los puntos $y_i$'s:
y si $A$ es de $rank$ completo (tiene $n+1$ columnas linealmente independientes) se calcula la factorización $QR$ de $A$ : $A = QR$ y entonces: $$A^TA\beta = A^Ty$$
y como $A=QR$ se tiene: $A^TA = (R^TQ^T)(QR)$ y $A^T = R^TQ^T$ por lo que:
y usando que $Q$ tiene columnas ortonormales:
Como $A$ tiene $n+1$ columnas linealmente independientes, la matriz $R$ es invertible por lo que $R^T$ también lo es y finalmente se tiene el sistema de ecuaciones lineales a resolver:
Caso regresión lineal
En el caso de la regresión lineal se ajusta un modelo de la forma: $f(x|\beta) = \beta_0 + \beta_1 x$ a los datos $(x_i,y_i)$'s $\forall i=0,1,\dots,m$.
Obs: En este caso se eligen $\phi_0(x) = 1$, $\phi_1(x) =x$. Y tenemos que estimar dos parámetros: $\beta_0, \beta_1$.
In [23]:
import matplotlib.pyplot as plt
import pprint
from scipy.linalg import solve_triangular
In [24]:
np.random.seed(1989) #para reproducibilidad
mpoints = 20
x = np.random.randn(mpoints)
y = -3*x + np.random.normal(2,1,mpoints)
In [25]:
plt.plot(x,y, 'r*')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Puntos ejemplo')
plt.show()
Con numpy podemos usar la función polyfit en el paquete de numpy para realizar el ajuste: (ver numpy.polyfit)
In [26]:
# el tercer argumento de polyfit especifica el grado del polinomio a ajustar.
#Usaremos ngrado = 1 pues queremos ajustar una recta
ngrado = 1
coeficientes = np.polyfit(x,y,ngrado)
In [27]:
#Una vez realizado el llamado a la función polyfit se regresan los coeficientes de x
#ordenados del mayor grado al menor.
pprint.pprint(coeficientes)
Entonces nuestro polinomio es: $$p_{1}(x) = -2.65x + 2.03$$
y así tenemos nuestras beta's ajustadas $\hat{\beta_0} = 2.03$, $\hat{\beta_1} = -2.65$.
Ahora nos gustaría graficar el modelo en el intervalo $[min(x),max(x)]$ con $min(x)$ la entrada con valor mínimo del numpy array $x$ y $max(x)$ su entrada con valor máximo.
Para lo anterior debemos obtener los valores ajustados al evaluar $p_1(x)$ los valores de $x$:
In [28]:
y_ajustadas_numpy = coeficientes[1] + coeficientes[0] * x
In [29]:
plt.plot(x, y_ajustadas_numpy, 'k-',x, y, 'r*')
plt.legend(['modelo lineal','datos'], loc='best')
plt.show()
In [30]:
#construimos a la matriz A:
A=np.ones((mpoints,2))
In [31]:
A[:,1] = x
In [32]:
A
Out[32]:
In [33]:
Q,R = np.linalg.qr(A)
In [34]:
#Resolvemos el sistema R*beta = Q^T*y
beta = solve_triangular(R,Q.T@y)
pprint.pprint(beta)
In [35]:
y_ajustadas_QR = A@beta
#obsérvese que la línea anterior es equivalente a realizar:
#y_ajustadas_QR = beta[0] + beta[1]*x
In [36]:
plt.plot(x, y_ajustadas_QR , 'k-',x, y, 'r*')
plt.legend(['modelo lineal','datos'], loc='best')
plt.show()
Usamos optimización numérica reescribiendo la función objetivo:
con $y \in \mathbb{R}^{20}, A \in \mathbb{R}^{20 \times 2}, \beta \in \mathbb{R}^{2 \times 1}$ y $A[i,:]$ $i$-ésimo renglón de $A$.
Planteamos el problema de optimización numérica:
$$\displaystyle \min_{\beta \in \mathbb{R}^n} ||y - A\beta||_2^2$$¿Solución?
Reescribimos $f_o$ como $f_o(\beta) = ||y-A\beta||_2^2= (y-A\beta)^T(y-A\beta) = y^Ty-2\beta^TA^Ty + \beta^TA^TA\beta$, observa que esta última expresión es un número $\mathbb{R}$.
Por lo anterior: $\nabla f_o(\beta) = -2A^Ty + 2A^TA\beta$. $\nabla f_o$ es un vector en $\mathbb{R}^{2 \times 1}$, de hecho sus entradas son $\nabla f_o(\beta) = \left [ \begin{array}{c} \frac{\partial f_o(\beta)}{\beta_0}\\ \frac{\partial f_o(\beta)}{\beta_1} \end{array} \right ] $ .
Si se plantea la ecuación no lineal: $\nabla f_o(\beta)= 0$ obtenemos una forma cerrada de la solución que está dada por $A^TA \beta=A^Ty$ (y son las ecuaciones normales!).
Consideramos métodos de descenso para resolver el problema de optimización.
La solución como la encontramos en el ejemplo del inicio, es considerar un método iterativo en el que tomamos un punto inicial $\beta^{(0)}$ y las actualizaciones se realizan:
y $t_{k}$ obtenida por búsqueda de línea con backtracking y es una cantidad positiva.
Obs: igualar la ecuación lineal $-2A^Ty + 2A^TA\beta$ a 0 es equivalente a resolver la ecuación lineal $-A^Ty + A^TA\beta=0$, por esto utilizamos esta última expresión que corresponde a una función objetivo $f_o = \frac{1}{2}y^Ty-\beta^TA^Ty + \frac{1}{2}\beta^TA^TA\beta$.
calculamos primero $-A^Ty$:
In [37]:
cte=-np.transpose(A)@y
In [38]:
cte
Out[38]:
el gradiente de $f_o$: $\nabla f_o(\beta) = -A^Ty + A^TA\beta$:
In [39]:
gf = lambda beta_fun: cte + np.transpose(A)@(A@beta_fun)
#observa que no hacemos la multiplicación (A^T*A)*beta, mejor hacemos
#primero A*beta y luego multiplicamos por A^T
Iniciamos con las iteraciones dadas por la fórmula de actualización tomando $\beta_0=(0,0)^T$:
In [40]:
beta_0 = np.array([0,0])
Comentario: más adelante veremos cómo obtener las $t_{k}$s de la búsqueda de línea por backtracking
In [41]:
t_0=.130
In [42]:
beta_1 = beta_0 - t_0*gf(beta_0)
In [43]:
beta_1
Out[43]:
In [44]:
t_1=.0625
In [45]:
beta_2 = beta_1 - t_1*gf(beta_1)
In [46]:
beta_2
Out[46]:
In [47]:
t_2 = .0625
In [48]:
beta_3 = beta_2 - t_2*gf(beta_2)
In [49]:
beta_3
Out[49]:
In [50]:
t_3 = .0625
In [51]:
beta_4 = beta_3 - t_3*gf(beta_3)
In [52]:
beta_4
Out[52]:
In [53]:
t_4 = .0625
In [54]:
beta_5 = beta_4 - t_4*gf(beta_4)
In [55]:
beta_5
Out[55]:
In [56]:
y_ajustadas_gradiente = A@beta_5
#obsérvese que la línea anterior es equivalente a realizar:
#y_ajustadas_gradiente = beta_5[0] + beta_5[1]*x
In [57]:
plt.plot(x, y_ajustadas_gradiente , 'k-',x, y, 'r*')
plt.legend(['modelo lineal','datos'], loc='best')
plt.show()
Otra opción es utilizar la información de segundo orden con la Hessiana y considerar una actualización:
para $k=0,1,2,\dots,$ (omitiendo la búsqueda de línea por backtracking). Recuérdese que no invertimos a la matriz y en su lugar resolvemos un sistema de ecuaciones lineales asociado. Para el ejemplo de mínimos cuadrados lineales se tiene que la Hessiana es:
Y el sistema de ecuaciones lineales asociado es: $\nabla ^2 f(\beta ^{(k)}) s_k = - \nabla f(\beta^{(k)})$:
In [58]:
s_0 = np.linalg.solve(np.transpose(A)@A,-gf(beta_0))
beta_1 = beta_0 + s_0
In [59]:
beta_1
Out[59]:
In [60]:
y_ajustadas_Newton = A@beta_1
#obsérvese que la línea anterior es equivalente a realizar:
#y_ajustadas_Newton = beta_1[0] + beta_1[1]*x
In [61]:
plt.plot(x, y_ajustadas_Newton , 'k-',x, y, 'r*')
plt.legend(['modelo lineal','datos'], loc='best')
plt.show()
Hasta este punto comparemos con lo obtenido previamente por QR o polyfit
Para descenso en gradiente con búsqueda de línea por backtracking se tiene un error relativo de:
In [62]:
np.linalg.norm(beta_5-beta)/np.linalg.norm(beta)
Out[62]:
Para descenso con dirección de Newton (sin búsqueda de línea por backtracking) se tiene un error relativo de::
In [63]:
np.linalg.norm(beta_1-beta)/np.linalg.norm(beta)
Out[63]:
Comentarios:
Algoritmo de descenso
Dado un punto inicial $x$ en $\text{dom}f$
Repetir el siguiente bloque para $k=0,1,2,...$
- Determinar una dirección de descenso $\Delta x$ (gradiente o Newton por ejemplo).
- Búsqueda de línea. Elegir un tamaño de paso $t > 0$.
- Hacer la actualización: $x = x + t\Delta x$.
hasta convergencia (satisfacer criterio de paro).
La Hessiana en mínimos cuadrados lineales como se vio en este ejemplo es $\nabla^2f(\beta^{(k)}) = A^TA$. No es recomendable construir la matriz $A^TA$ por dos razones principalmente:
Si $A$ es grande (más de $10^8$ entradas) hacer la operación $A^TA$ costará gran cantidad de tiempo por el número de operaciones involucradas (alrededor de $\mathcal{O}(mn^2)$).
Si $A$ tiene un número de condición $c>0$ se prueba que $A^TA$ tiene un número de condición $c^2$ bajo la norma $2$ . Esto como se revisó en la nota 1.3.Condicion_de_un_problema_y_estabilidad_de_un_algoritmo afecta en la precisión y exactitud al resolver el problema de encontrar la solución al sistema de ecuaciones lineales $\nabla ^2 f(\beta^{(k)})s_k = - \nabla f(\beta^{(k)})$.
Por lo anterior se sugiere resolver el sistema de ecuaciones lineales $\nabla ^2 f(\beta ^{(k)}) s_k = - \nabla f(\beta^{(k)})$ dado por: $A^TA s_k = -(-A^Ty + A^TA\beta) $ que resulta de los mínimos cuadrados con la factorización $QR$:
donde se asume de la segunda igualdad a la tercera igualdad que $A$ es de rank completo y por lo tanto $R$ y $R^T$ son no singulares (son invertibles).
In [64]:
Q,R = np.linalg.qr(A)
In [65]:
s_0 = np.linalg.solve(R,np.transpose(Q)@y-R@beta_0)
beta_1 = beta_0 + s_0
In [66]:
beta_1
Out[66]:
In [67]:
y_ajustadas_Newton = A@beta_1
#obsérvese que la línea anterior es equivalente a realizar:
#y_ajustadas_Newton = beta_1[0] + beta_1[1]*x
In [68]:
plt.plot(x, y_ajustadas_Newton , 'k-',x, y, 'r*')
plt.legend(['modelo lineal','datos'], loc='best')
plt.show()
Obsérvese que esta forma es igual a la vista por factorización QR anteriormente ya que se resuelve el sistema de ecuaciones normales $A^TA \beta = A^Ty$.
Para la formulación del método consideramos a la función objetivo $f_o$ escrita en la forma:
$$f_o(\beta) = \frac{1}{2}\displaystyle \sum_{i=1}^{20} (y_i - A[i,:]^T\beta)^2$$con $A[i,:]$ $i$-ésimo renglón de $A$ y obsérvese que dividimos por $2$ como antes. Calculamos el gradiente de $f_o$:
$$\nabla f_o (\beta) = -\displaystyle \sum_{i=1}^{20} (y_i - A[i,:]^T\beta)A[i,:]$$obs: La expresión anterior considera $A[i,:] \in \mathbb{R}^{2 \times 1}$, $A[i,:]=\left[ \begin{array}{c} 1\\ x_i \end{array} \right ] $ $\forall i=1,\dots, 20$ con $x_i$ dato dado en el ejemplo de mínimos cuadrados lineales.
En el método estocástico utilizando el gradiente (nombrado descenso en gradiente estocástico), utilizamos un punto inicial $\beta^{(0)}$ y en lugar de usar la actualización:
$$ \begin{eqnarray} \beta^{(k+1)}&=&\beta^{(k)} - t_{k}\nabla f\left(\beta^{(k)}\right) \nonumber\\ &=& \beta^{(k)} - t_{k}\left(-\displaystyle \sum_{i=1}^{20} (y_i - A[i,:]^T\beta^{(k)})A[i,:]\right)\nonumber\\ \end{eqnarray} $$para $k=0,1,2,\dots$ y $t_{k}$ cantidad positiva, obtenida por búsqueda de línea con backtracking. Las actualizaciones en descenso en gradiente estocástico son:
para $k=0,1,2,\dots$ con $\nabla f_{i_{k}}(\cdot)$ el gradiente calculado a partir de una muestra extraída del conjunto de índices de los renglones: $\{1,2,\dots,m\}$.
Comentarios:
La cantidad $\eta_{k}$ en este contexto se le nombra tasa de aprendizaje y al igual que en el método de descenso en gradiente es una cantidad positiva.
La muestra de índices para calcular $\nabla f_{i_{k}}(\cdot)$ puede contener uno o más índices. Si contiene un único índice por ejemplo $i_{k} = \{i\}$ con $i=1,\cdots, m$ se tiene una actualización de la forma:
Para el ejemplo numérico de mínimos cuadrados lineales, iniciamos con las iteraciones dadas por la fórmula de actualización tomando $\beta_0=(0,0)^T$ y $\eta_{k} = 0.1 \forall k=0,1,2,\dots$:
In [69]:
beta_0 = np.array([0,0])
eta_0 = .1
Elegimos un índice de los renglones de $A$ de forma pseudoaleatoria y hacemos la actualización:
In [70]:
np.random.seed(1989) #para reproducibilidad
i = np.random.randint(mpoints, size=1)
In [71]:
beta_1 = beta_0 + eta_0*(y[i]-A[i,:].dot(beta_0))*A[i,:]
In [72]:
beta_1
Out[72]:
El algoritmo es:
Dado un punto inicial $\beta_0$.
Repetir el siguiente bloque para $k=0,1,2,\dots$ :
- Extraer una muestra aleatoria $i_{k}$ del conjunto de índices de los renglones de $A$: $\{1,2,\dots,m\}$.
- Calcular la tasa de aprendizaje $\eta_{k}$ (positiva).
- Calcular la actualización $\beta^{(k+1)} = \beta^{(k)} + \eta_{k}\nabla f_{i_{k}}\left(\beta^{(k)}\right)$
hasta convergencia.
Midamos el error relativo con este algoritmo utilizando sólo un índice $i$ extraído al azar del conjunto de índices de renglones de $A$ para el ejemplo de mínimos cuadrados lineales: $\beta^{(k+1)} = \beta^{(k)} + \eta_{k}(y_i - A[i,:]^T\beta^{(k)})A[i,:]$
In [73]:
def error_relativo(aprox, obj):
return np.linalg.norm(aprox-obj)/np.linalg.norm(obj)
Definamos puntos iniciales
In [76]:
maxiter=5
tol=1e-2
err=error_relativo(beta_0,beta)
beta_k = beta_0
k=1
eta_k = .1 #constante para todas las iteraciones
In [77]:
err
Out[77]:
In [78]:
while(err>tol and k <= maxiter):
i = int(np.random.randint(mpoints, size=1))
beta_k = beta_k+eta_k*(y[i]-A[i,:].dot(beta_k))*A[i,:]
k+=1
eta_k=.1
err=error_relativo(beta_k,beta)
In [79]:
err
Out[79]:
tenemos un error relativo de $92\%$ con $5$ iteraciones.
In [80]:
beta_k
Out[80]:
Consideremos una muestra más grande de índices para calcular $\nabla f_{i_{k}}(\cdot)$ de modo que las actualizaciones ahora son de la forma:
$$ \begin{eqnarray} \beta^{(k+1)} &=& \beta^{(k)} - \eta_{{k}}\nabla f_{i_{k}}\left(\beta^{(k)}\right) \nonumber\\ &=& \beta^{(k)} - \eta_{k}\displaystyle \sum_{i \in i_{k}} (-(y_i - A[i,:]^T\beta^{(k)})A[i,:]) \nonumber\\ &=& \beta^{(k)} + \eta_{k}\displaystyle \sum_{i \in i_{k}}(y_i - A[i,:]^T\beta^{(k)})A[i,:] \nonumber\\ \end{eqnarray} $$para $k=0,1,2,\dots$ e $i_{k}$ una muestra de índices extraída de forma aleatoria de los índices de los renglones de $A$.
Usamos $4$ iteraciones con una muestra de índices de tamaño $5$. Puntos iniciales:
In [81]:
maxiter=4
tol=1e-2
err=error_relativo(beta_0,beta)
beta_k = beta_0
k=1
eta_k = .1 #constante para todas las iteraciones
m_sample = 5
In [82]:
err
Out[82]:
In [83]:
while(err>tol and k <= maxiter):
idx = np.random.choice(mpoints, m_sample,replace=False) #muestra de tamaño m_sample sin reemplazo
beta_k = beta_k+eta_k*(y[idx]-A[idx,:]@beta_0)@A[idx,:]
k+=1
eta_k=.1
err=error_relativo(beta_k,beta)
In [84]:
err
Out[84]:
Obsérvese que tenemos un error relativo de $51\%$ con $4$ iteraciones, una muestra de tamaño $5$ y es menor al que se obtuvo con una muestra de tamaño 1 de índices.
In [85]:
beta_k
Out[85]:
Comentarios: Algunos comentarios generales para el enfoque por batch o lote y el estocástico que podemos realizar son:
El enfoque por batch:
El enfoque estocástico:
*Una razón intuitiva del por qué el método estocástico de descenso en gradiente podría ganarle en velocidad a una formulación en batch es la siguiente: supóngase que un conjunto de entrenamiento $\mathcal{S}$ consistiera de renglones redundantes. Imaginemos que $\mathcal{S}$ consta de $10$ repeticiones de un conjunto más pequeño $\mathcal{S}_\text{sub}$. Entonces ejecutar un método batch sería $10$ veces más costoso que si solamente tuviéramos una única copia de $\mathcal{S}_\text{sub}$. El método estocástico de descenso en gradiente tendría el mismo costo en cualquiera de los dos escenarios ($10$ copias vs $1$ copia) pues elige con misma probabilidad renglones del conjunto $\mathcal{S}_\text{sub}$. En la realidad un conjunto de entrenamiento no contiene duplicados exactos de la información pero en muchas aplicaciones large-scale existe (de forma aproximada) información redundante.
Aplicaciones de machine learning conducen al planteamiento de problemas de optimización convexa y no convexa.
Ejemplo de problemas convexos los encontramos en la aplicación de clasificación de textos en donde se desea asignar un texto a clases definidas de acuerdo a su contenido (p.ej. determinar si un docu de texto es sobre política). La aplicación anterior puede formularse utilizando funciones de pérdida convexas.
Como ejemplos de aplicaciones en el ámbito de la optimización no convexa están el reconocimiento de voz y reconocimiento de imágenes. El uso de redes neuronales profundas* ha tenido muy buen desempeño en tales aplicaciones haciendo uso de cómputo en la GPU, ver 2.3.CUDA, ImageNet Classification with Deep Convolutional Neural Networks, 2012: A Breakthrough Year for Deep Learning. En este caso se utilizan funciones objetivo no lineales y no convexas.
*Los tipos de redes neuronales profundas, deep neural networks, que han sido mayormente usadas en el inicio del siglo XXI son las mismas que las que eran populares en los años $90$'s. El éxito de éstos tipos y su uso primordialmente se debe a la disponibilidad de larger datasets y mayores recursos computacionales.
Desde los $40$'s se han desarrollado algoritmos para resolver problemas de optimización, se han analizado sus propiedades y se han desarrollado buenas implementaciones de software. Sin embargo, una clase de problemas de optimización en los que encontramos métodos efectivos son los convexos.
En el módulo IV del curso nos enfocamos mayormente a métodos numéricos para resolver problemas convexos principalmente por lo anterior y porque métodos para optimización no convexa utilizan parte de la teoría de convexidad desarrollada en optimización convexa. Además un buen número de problemas de aprendizaje utilizan funciones de pérdida convexas.
En lo que continúa se considera $f_0 = f_o$ (el subíndice "0" y el subíndice "o" son iguales)
con $f_i: \mathbb{R}^n \rightarrow \mathbb{R}$ $\forall i=0,\dots,m$, $h_i: \mathbb{R}^n \rightarrow \mathbb{R}$, $\forall i=1,\dots,p$. $f_i$ son las restricciones de desigualdad, $h_i$ son las restricciones de igualdad.
(en github aparecen los símbolos de intersección más abajo de lo que deberían aparecer...abriendo con jupyter el notebook esto no pasa...)
Comentarios:
Un punto $x \in \mathcal{D}$ se nombra factible si satisface las restricciones de igualdad y desigualdad. El conjunto de puntos factibles se nombra conjunto de factibilidad.
El problema anterior se nombra problema de optimización factible si existe al menos un punto factible, si no se cumple lo anterior entonces es infactible.
En lo que continúa se asumen todos los puntos en el dominio $\mathcal{D}$.
El valor óptimo del problema se denota como $p^*$. En notación matemática es:
$$p^* = \inf\{f_o(x) | f_i(x) \leq 0, \forall i=1,\dots,m, h_i(x) = 0 \forall i=1,\dots,p\}$$Comentarios:
Si el problema es infactible entonces $p^* = \infty$.
Si $\exists x_k$ factible tal que $f_o(x_k) \rightarrow -\infty$ para $k \rightarrow \infty$ entonces $p^*=-\infty$ y se nombra problema de optimización no acotado por debajo.
$x^*$ es punto óptimo si es factible y $f_o(x^*) = p^*$. El conjunto de óptimos se nombra conjunto óptimo y se denota:
$$X_{\text{opt}} = \{x | f_i(x) \leq 0 \forall i=1,\dots,m, h_i(x) =0 \forall i=1,\dots,p, f_o(x) = p^*\}$$Comentarios:
La propiedad de un punto óptimo $x^*$ es que si $z$ satisface las restricciones $f_i(z) \leq 0$ $\forall i=1,...,m$, $h_i(z)=0$ $\forall i=1,..,p$ se tiene: $f(x^*) \leq f(z)$. Es óptimo estricto si $z$ satisface las restricciones y $f_o(x^*) < f_o(z)$.
Si existe un punto óptimo se dice que el valor óptimo se alcanza y por tanto el problema de optimización tiene solución, es solvable.
Si $X_{\text{opt}} = \emptyset$ se dice que el valor óptimo no se alcanza. Obsérvese que para problemas no acotados nunca se alcanza el valor óptimo.
Si $x$ es factible y $f_o(x) \leq p^* + \epsilon$ con $\epsilon >0$, $x$ se nombra $\epsilon$-subóptimo y el conjunto de puntos $\epsilon$-subóptimos se nombra conjunto $\epsilon$-subóptimo.
Un punto factible $x^*$ se nombra óptimo local si $\exists R > 0$ tal que:
$$f_o(x^*) = \inf \{f_o(z) | f_i(z) \leq 0 \forall i=1,\dots,m, h_i(z) = 0 \forall i=1,\dots, p, ||z-x||_2 \leq R\}.$$Así, $x^*$ resuelve:
$$\displaystyle \min f_o(z)$$Obs: la palabra óptimo se utiliza para óptimo global, esto es, no consideramos la última restricción $||z-x||_2 \leq R$ en el problema de optimización y exploramos en todo el $\text{dom}f$.
Si $x$ es factible y $f_i(x)=0$ entonces la restricción de desigualdad $f_i(x) \leq 0$ se nombra restricción activa en $x$. Se nombra inactiva en $x$ si $f_i(x) <0$ para alguna $i$.
Comentario: Las restricciones de igualdad, $h_i(x)$, siempre son activas en el conjunto factible.
Una restricción se nombra restricción redundante si al quitarla el conjunto factible no se modifica.
donde: $f_i$ son convexas $\forall i=0,2,\dots,m$ y $h_i$ es afín* $\forall i =1,\dots,p$.
*Una función afín es de la forma $h(x) = Ax+b$ con $A \in \mathbb{R}^{p \times n}$ y $b \in \mathbb{R}^p$. En la definición anterior $h_i(x) = a_i^Tx-b_i$ con $a_i \in \mathbb{R}^n$, $b_i \in \mathbb{R}$ $\forall i=1,\dots,p$ y geométricamente $h_i(x)$ es un hiperplano en $\mathbb{R}^n$.
Comentarios: Lo siguiente se puede verificar con propiedades teóricas:
*Un conjunto $\alpha$-subnivel es de la forma $\{x \in \text{dom}f | f(x) \leq \alpha\}$. Un conjunto subnivel contiene las curvas de nivel de $f$, ver Level set:
Si $f$ es convexa el conjunto subnivel es un conjunto convexo.
El conjunto óptimo y los conjuntos $\epsilon$-subóptimos son convexos.
Si la función objetivo $f_o$ es fuertemente convexa entonces el conjunto óptimo contiene a lo más un punto.
Si en el problema anterior se tiene que maximizar una $f_o$ función objetivo cóncava y se tienen misma forma estándar del problema anterior y $f_i$ convexa, $h_i$ afín entonces también se nombra al problema como problema de optimización convexa. Todos los resultados, conclusiones y algoritmos desarrollados para los problemas de minimización son aplicables para maximización. En este caso se puede resolver un problema de maximización al minimizar la función objetivo $-f_o$ que es convexa.
Sea $f:\mathbb{R}^n \rightarrow \mathbb{R}$ una función con el conjunto $\text{dom}f$ convexo. $f$ se nombra convexa (en su $\text{dom}f$) si $\forall x,y \in \text{dom}f$ y $\theta \in [0,1]$ se cumple:
$$f(\theta x + (1-\theta) y) \leq \theta f(x) + (1-\theta)f(y).$$Comentarios:
$\text{dom}f$ es convexo $\therefore$ $\theta x + (1-\theta)y \in \text{dom}f$
La convexidad de $f$ se define para $\text{dom}f$ aunque para casos en particular se detalla el conjunto en el que $f$ es convexa.
Si la desigualdad se cumple de forma estricta $\forall x \neq y$ $f$ se nombra estrictamente convexa.
$f$ es cóncava si $-f$ es convexa y estrictamente cóncava si $-f$ es estrictamente convexa. Otra forma de definir concavidad es con una desigualdad del tipo:
y mismas definiciones para $x,y, \theta$ que en la definición de convexidad.
Sean $x_1, x_2 \in \mathbb{R}^n$ con $x_1 \neq x_2$. Entonces el punto:
$$y = \theta x_1 + (1-\theta)x_2$$con $\theta \in \mathbb{R}$ se encuentra en la línea que pasa por $x_1$ y $x_2$. $\theta$ se le nombra parámetro y si $\theta \in [0,1]$ tenemos un segmento de línea:
Comentarios:
$y = \theta x_1 + (1-\theta)x_2 = x_2 + \theta(x_1 -x_2)$ y esta última igualdad se interpreta como "$y$ es la suma del punto base $x_2$ y la dirección $x_1-x_2$ escalada por $\theta$".
Si $\theta=0$ entonces $y=x_2$. Si $\theta \in [0,1]$ entonces $y$ se "mueve" en la dirección $x_1-x_2$ hacia $x_1$ y si $\theta>1$ entonces $y$ se encuentra en la línea "más allá" de $x_1$:
El punto enmedio entre $x_1$ y $x_2$ tiene $\theta=\frac{1}{2}$.
Un conjunto $\mathcal{C}$ es convexo si el segmento de línea entre cualquier par de puntos de $\mathcal{C}$ está completamente contenida en $\mathcal{C}$. Esto se escribe matemáticamente como:
$$\theta x_1 + (1-\theta) x_2 \in \mathcal{C} \forall \theta \in [0,1], \forall x_1, x_2 \in \mathcal{C}.$$Ejemplos gráficos de conjuntos convexos:
Ejemplos gráficos de conjuntos no convexos:
Comentarios:
El punto $\displaystyle \sum_{i=1}^k \theta_i x_i$ con $\displaystyle \sum_{i=1}^k \theta_i=1$, $\theta_i \geq 0 \forall i=1,\dots,k$ se nombra combinación convexa de los puntos $x_1, x_2, \dots, x_k$. Una combinación convexa de los puntos $x_1, \dots, x_k$ puede pensarse como una mezcla o promedio ponderado de los puntos, con $\theta_i$ la fracción $\theta_i$ de $x_i$ en la mezcla.
Un conjunto es convexo si y sólo si contiene cualquier combinación convexa de sus puntos.
Obs: por tanto las funciones lineales también son convexas y cóncavas.
*Una matriz $A$ es positiva semidefinida si $x^TAx \geq 0$ $\forall x \in \mathbb{R}^n - \{0\}$. Si se cumple de forma estricta la desigualdad anterior entonces $A$ es positiva definida. Con los eigenvalores podemos caracterizar a las matrices definidas y semidefinidas positivas: $A$ es semidefinida positiva si y sólo si los eigenvalores de $T=\frac{A+A^T}{2}$ son no negativos. Es definida positiva si y sólo si los eigenvalores de $T$ son positivos. Los conjuntos de matrices que se utilizan para definir a matrices semidefinidas positivas y definidas positivas son $\mathbb{S}_{+}^n$ y $\mathbb{S}_{++}^n$ respectivamente.
Obs: $f$ es estrictamente convexa si y sólo si $P \in \mathbb{S}_{++}^n$. $f$ es cóncava si y sólo si $P \in -\mathbb{S}_+^n$.
Exponenciales: $f: \mathbb{R} \rightarrow \mathbb{R}$, $f(x) = e^{ax}$ para cualquier $a \in \mathbb{R}$ es convexa en su dominio: $\mathbb{R}$.
Potencias: $f: \mathbb{R} \rightarrow \mathbb{R}$, $f(x)=x^a$:
Potencias del valor absoluto: $f: \mathbb{R} \rightarrow \mathbb{R}$, $f(x)=|x|^p$ con $p \geq 1$ es convexa en $\mathbb{R}$.
Logaritmo: $f: \mathbb{R} \rightarrow \mathbb{R}$, $f(x) = \log(x)$ es cóncava en su dominio: $\mathbb{R}_{++}$.
Entropía negativa: $f(x) = \begin{cases} x\log(x) &\text{ si } x > 0 ,\\ 0 &\text{ si } x = 0 \end{cases}$ es estrictamente convexa en su dominio: $\mathbb{R}_+$.
Normas: cualquier norma es convexa en su dominio.
Función máximo: $f: \mathbb{R}^{n} \rightarrow \mathbb{R}$, $f(x) = \max\{x_1,\dots,x_n\}$ es convexa.
Función log-sum-exp: $f: \mathbb{R}^{n} \rightarrow \mathbb{R}$, $f(x)=\log\left(\displaystyle \sum_{i=1}^ne^{x_i}\right)$ es convexa en su dominio: $\mathbb{R}^n$.
La media geométrica: $f: \mathbb{R}^{n} \rightarrow \mathbb{R}$, $f(x) = \left(\displaystyle \prod_{i=1}^n x_i \right)^\frac{1}{n}$ es cóncava en su dominio: $\mathbb{R}_{++}^n$.
Función log-determinante: $f: \mathbb{S}^{n} \rightarrow \mathbb{R}^n$, $f(x) = \log(\det(X))$ es cóncava en su dominio: $\mathbb{S}_{++}^n$.
Se sugiere revisar 1.4.Polinomios_de_Taylor_y_diferenciacion_numerica y 1.3.Condicion_de_un_problema_y_estabilidad_de_un_algoritmo como recordatorio de definiciones. En particular las definiciones de primera y segunda derivada, gradiente y Hessiana para la primer nota y la definición de número de condición de una matriz para la segunda.
Sobre funciones convexas/cóncavas
Si se cumple de forma estricta la desigualdad $f$ se nombra estrictamente convexa. También si su $\text{dom}f$ es convexo y se tiene la desigualdad en la otra dirección "$\leq$" entonces $f$ es cóncava.
Geométricamente este resultado se ve como sigue para $\nabla f(x) \neq 0$:
y el hiperplano $f(x) + \nabla f(x)^T(y-x)$ se nombra hiperplano de soporte para la función $f$ en el punto $(x,f(x))$. Obsérvese que si $\nabla f(x)=0$ se tiene $f(y) \geq f(x) \forall y \in \text{dom}f$ y por lo tanto $x$ es un mínimo global de $f$.
*Para una variable: $f: \mathbb{R} \rightarrow \mathbb{R}$, la hipótesis del enunciado anterior es que la segunda derivada sea positiva. El recíproco no es verdadero, para ver esto considérese $f(x)=x^4$ la cual es estrictamente convexa en $\text{dom}f$ pero su segunda derivada en $0$ no es positiva.
Sobre problemas de optimización
Si $f$ es diferenciable y $x^*$ es óptimo entonces $\nabla f(x^*) = 0$.
Si $f \in \mathcal{C}^2(\text{domf})$ y $x^*$ es mínimo local entonces $\nabla^2 f(x^*)$ es una matriz simétrica semidefinida positiva.
Si $f \in \mathcal{C}^2(\text{domf})$, $\nabla f(x^*)=0$ y $\nabla^2f(x^*)$ es una matriz simétrica definida positiva entonces $x^*$ es mínimo local estricto.
Una propiedad fundamental de un óptimo local en un problema de optimización convexa es que también es un óptimo global.
Si $f_o$ en un problema de optimización convexa es diferenciable y $X$ es el conjunto de factibilidad entonces $x$ es óptimo si y sólo si $x \in X$ y $\nabla f_o(x)^T(y-x) \geq 0$ $\forall y \in X$. Si en lo anterior se considera como conjunto de factibilidad $X = \text{dom}f_o$ (que es un problema sin restricciones) la propiedad se reduce a la condición: $x$ es óptimo si y sólo si $\nabla f_o(x) = 0$.
Geométricamente el resultado anterior se visualiza para $\nabla f(x) \neq 0$ y $-\nabla f(x)$ apuntando hacia la dirección dibujada:
Comentario: Por los resultados anteriores los métodos de optimización que revisaremos en el módulo IV buscarán resolver la ecuación no lineal $\nabla f_o(x)=0$. Dependiendo del número de soluciones de la ecuación $\nabla f_o(x)=0$ se tienen situaciones distintas. Por ejemplo, si no tiene solución entonces el/los óptimos no se alcanza(n) pues el problema puede no ser acotado por debajo o si existe el óptimo éste puede no alcanzarse. Por otro lado, si la ecuación tiene múltiples soluciones entonces cada solución es un mínimo de $f_o$.
Una función $f:\mathbb{R}^n \rightarrow \mathbb{R}$ tal que $f \in \mathcal{C}^2(\text{dom}f)$ se nombra fuertemente convexa en el conjunto convexo $\mathbb{S} \neq \emptyset$ si existe $m>0$ tal que $\nabla^2 f(x) - mI$ es simétrica semidefinida positiva $\forall x \in \mathbb{S}$.
Comentario: si una función es fuertemente convexa se puede probar que:
$f(y) \geq f(x) + \nabla f(x)^T(y-x) + \frac{m}{2}||y-x||_2^2 \forall x,y \in \mathbb{S}$. Por esto si $f$ es fuertemente convexa en $\mathbb{S}$ entonces es estrictamente convexa en $\mathbb{S}$. También esta desigualdad indica que la diferencia entre la función de $y$, $f(y)$, y la función lineal en $y$ $f(x) + \nabla f(x)^T(y-x)$ (Taylor a primer orden) está acotada por debajo por una cantidad cuadrática.
El número de condición bajo la norma 2 de $\nabla ^2 f$ está acotado por arriba por el cociente $\frac{\lambda_\text{max}(\nabla^2 f(x))}{\lambda_\text{min}(\nabla^2 f(x))} \forall x \in \mathbb{S}$.
Si $f,g:\mathbb{R}^n \rightarrow \mathbb{R}$ con $f,g \in \mathcal{C}^2$ respectivamente en sus dominios y $\alpha_1, \alpha_2 \in \mathbb{R}$, $A \in \mathbb{R}^{n \times n}$, $b \in \mathbb{R}^n$ son fijas. Diferenciando con respecto a la variable $x \in \mathbb{R}^n$ se tiene:
| linealidad | $\nabla(\alpha_1 f(x) + \alpha_2 g(x)) = \alpha_1 \nabla f(x) + \alpha_2 \nabla g(x)$ |
| producto | $\nabla(f(x)g(x)) = \nabla f(x) g(x) + f(x) \nabla g(x)$ |
| producto punto | $\nabla(b^Tx) = b$ |
| cuadrático | $\nabla(x^TAx) = (A+A^T)x$ |
| segunda derivada | $\nabla^2(Ax)=A$ |
Referencias:
S. P. Boyd, L. Vandenberghe, Convex Optimization, Cambridge University Press, 2009.
L. Bottou, F. E. Curtis, J Nocedal, Optimization Methods for Large-Scale Machine Learning, SIAM, 2018.