import numpy as np
import sympy
from sympy.tensor.array import derive_by_array
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
from scipy.optimize import fmin
import pandas as pd
import cvxpy as cp
from pytest import approx

np.set_printoptions(precision=3, suppress=True)

"""
Explanation: (ADBLUCO)=
3.2 Algoritmos de descenso y búsqueda de línea en Unconstrained Convex Optimization (UCO)
```{admonition} Notas para contenedor de docker:
Comando de docker para ejecución de la nota de forma local:
nota: cambiar &lt;ruta a mi directorio&gt; por la ruta de directorio que se desea mapear a /datos dentro del contenedor de docker y &lt;versión imagen de docker&gt; por la versión más actualizada que se presenta en la documentación.
docker run --rm -v &lt;ruta a mi directorio&gt;:/datos --name jupyterlab_optimizacion -p 8888:8888 -d palmoreck/jupyterlab_optimizacion:&lt;versión imagen de docker&gt;
password para jupyterlab: qwerty
Detener el contenedor de docker:
docker stop jupyterlab_optimizacion
Documentación de la imagen de docker palmoreck/jupyterlab_optimizacion:&lt;versión imagen de docker&gt; en liga.
```

Nota generada a partir de liga.
```{admonition} Al final de esta nota la comunidad lectora:
:class: tip


Comprenderá el uso de la información de primer y segundo orden para resolver problemas small scale de minimización de funciones convexas mediante los algoritmos general de descenso y de búsqueda de línea por backtracking.


Aprenderá la importancia y relación con ecuaciones no lineales al resolver los problemas que involucran aproximar mínimos locales de funciones.


```
En esta nota consideramos problemas de optimización small scale sin restricciones. Aunque el término small scale es ambiguo pues depende la máquina en la que se realice el cómputo e involucra el número de variables o parámetros y cantidad de almacenamiento para datos, tomamos como small scale aquel problema de optimización sin restricciones en el que se tiene un número de variables del orden menor o igual a $10^3$.
Los problemas de optimización sin restricciones, Unconstrained Optimization (UO), son de la forma:
$$\min f_o(x)$$
con $f_o:\mathbb{R}^n \rightarrow \mathbb{R}$ función objetivo.
```{admonition} Comentario
Si la función objetivo es convexa se le nombra problema de optimización convexa sin restricciones,Unconstrained Convex Optimization (UCO).
```
Ejemplos de problemas de optimización small scale
En optimización la búsqueda del (o los) óptimo(s) involucran el cálculo de información de primer o segundo orden, ver {ref}Definición de función, continuidad y derivada &lt;FCD&gt;, de la función $f_o$ de acuerdo a lo revisado en los {ref}resultados útiles&lt;RESUT&gt;. Tal información para problemas small scale es calculada utilizando todos los datos en un enfoque por batch o lote.
Ejemplo
$$\displaystyle \min_{x \in \mathbb{R}^2} x_1^4+2x_1^2x_2+x_2^2$$
End of explanation
"""


x1, x2 = sympy.symbols("x1, x2")

fo_sympy = x1**4+2*x1**2*x2+x2**2

sympy.pprint(fo_sympy)

gf = derive_by_array(fo_sympy, (x1, x2))

sympy.pprint(gf)

"""
Explanation: Los candidatos a ser mínimos los encontramos al calcular el gradiente de $f_o$ que podemos calcular con SymPy
End of explanation
"""


def fo_numpy(x):
    return x[0]**4 + 2*x[0]**2*x[1] + x[1]**2

x1_plot,x2_plot = np.meshgrid(np.linspace(-2,2,100), np.linspace(-2,2,100))
z_plot = x1_plot**4 + 2*x1_plot**2*x2_plot + x2_plot**2

x1_np = 0
x2_np = 0 
z_np = fo_numpy([x1_np, x2_np])
point = (x1_np, x2_np, z_np)

print(point)

# Create the figure
fig = plt.figure()
# Add an axes
ax = fig.gca(projection='3d')
ax.plot_surface(x1_plot, x2_plot, z_plot, alpha=0.2)
ax.scatter(point[0], point[1], point[2], color='green')
plt.title("$f_o(x) = x_1^4+2x_1^2x_2+x_2^2$")
plt.show()

"""
Explanation: y plantear:
$$
\nabla f_o(x) = 
\left [
\begin{array}{c}
4x_1^3+4x_1x_2\
2x_1^2+2x_2
\end{array}
\right ]=0
$$
la cual es una ecuación de dos variables y dos incógnitas no lineal. Resolviendo para $x_2$ se obtiene la relación: $x_2 = -x_1^2$. Entonces todos los puntos con coordenadas $x = (x_1, x_2)$ que satisfacen tal relación cumplen $\nabla f_o(x) = 0$. ¿Todos serán mínimos locales?
Gráfica de la superficie $f_o$
End of explanation
"""


x1_plot,x2_plot = np.meshgrid(np.linspace(-2,2,100), np.linspace(-4, 1,100))
z_plot = x1_plot**4 + 2*x1_plot**2*x2_plot + x2_plot**2

plt.contour(x1_plot,x2_plot,z_plot)
plt.scatter(point[0], point[1], color="green")
plt.title("Curvas de nivel de $f_o$")
plt.show()

"""
Explanation: Gráfica de las curvas de nivel de $f_o$
End of explanation
"""


def eqn(x):
    x1,x2=x
    return [4*x1**3+4*x1*x2, 2*x1**2+2*x2]

"""
Explanation: Resolvamos con scipy.optimize.fsolve la ecuación no lineal $\nabla f_o(x) = 0$
End of explanation
"""


root1 = fsolve(eqn, (1, 1))

root2 = fsolve(eqn, (-1, 1))

root3 = fsolve(eqn, (2, 0))

dic_roots = {"root1": root1,
             "root2": root2,
             "root3": root3}

"""
Explanation: ```{margin}
Elegimos diferentes puntos iniciales.
```
End of explanation
"""


print(pd.DataFrame(dic_roots))

"""
Explanation: ```{margin}
Obsérvese que los puntos root1, root2, root3 satisfacen la relación $x_2 = -x_1^2$. 
```
End of explanation
"""


gf_eval = lambda x: np.array([partial_derivative.subs({"x1": x[0], 
                                                       "x2": x[1]}) for partial_derivative in gf],
                             dtype=float)

dic = {"root1": gf_eval(root1),
       "root2": gf_eval(root2),
       "root3": gf_eval(root3)}

"""
Explanation: Al evaluar el gradiente en cada punto obtenemos cero (o cercano a cero):
End of explanation
"""


print(pd.DataFrame(dic).round(3))

"""
Explanation: ```{margin}
Los puntos root1, root2 y root3 resuelven la ecuación no lineal $\nabla f(x) = 0$ .
```
End of explanation
"""


Hf = derive_by_array(gf, (x1, x2))

sympy.pprint(Hf)

"""
Explanation: ¿Cómo podemos identificar si son mínimos? ... usamos la Hessiana de $f_o$
End of explanation
"""


Hf_eval = lambda x: np.array([second_partial_derivative.subs({"x1": x[0], 
                                                              "x2": x[1]}) for second_partial_derivative in Hf],
                             dtype=float)

Hf_root1 = Hf_eval(root1)
Hf_root2 = Hf_eval(root2)
Hf_root3 = Hf_eval(root3)

"""
Explanation: y revisamos eigenvalores de la Hessiana evaluada en los puntos root1, root2, root3
End of explanation
"""


print(np.linalg.eigvals(Hf_root1))

"""
Explanation: ```{margin}
La Hessiana en root1 es semidefinida positiva por lo que no podemos concluir que sea mínimo local de $f_o$.
```
End of explanation
"""


print(np.linalg.eigvals(Hf_root2))

"""
Explanation: ```{margin}
La Hessiana en root2 es semidefinida positiva por lo que no podemos concluir que sea mínimo local de $f_o$.
```
End of explanation
"""


print(np.linalg.eigvals(Hf_root3))

"""
Explanation: ```{margin}
La Hessiana en root3 es semidefinida positiva por lo que no podemos concluir que sea mínimo local de $f_o$.
```
End of explanation
"""


sympy.pprint(sympy.Matrix(Hf).det())

"""
Explanation: Tenemos un criterio para $2$ dimensiones:
```{admonition} Comentario
Sea $f \in \mathcal{C}^2(\text{intdom}f)$, $\det(\nabla^2 f(x))$ determinante de la Hessiana y $x \in \mathbb{R}^2$ un punto crítico o estacionario de $f$:


Si $\frac{\partial^2f(x)}{\partial x_1^2} < 0$ y $\det(\nabla^2 f(x)) >0$ entonces $f$ tiene un máximo local en $x$.


Si $\frac{\partial^2f(x)}{\partial x_1^2} > 0$ y $\det(\nabla^2 f(x)) >0$ entonces $f$ tiene un mínimo local en $x$.


Si $\det(\nabla^2 f(x)) < 0$ entonces $f$ tiene un punto silla o saddle point en $x$.


Si $\det(\nabla^2 f(x)) = 0$ no podemos concluir si $x$ es extremo.


```
```{admonition} Observaciones
:class: tip


Al determinante de la Hessiana de $f$ se le nombra Hessiano de $f$.


Lo anterior es un caso particular de los resultados descritos en {ref}sobre puntos críticos &lt;SPCRITICOS&gt;.


```
En el ejemplo el Hessiano es:
End of explanation
"""


res_fmin = fmin(fo_numpy, [1,1])

"""
Explanation: ```{margin}
Los tres puntos root1, root2 y root3 satisfacen $x_2 = -x_1^2$.
```
El cual se anula justo en los puntos que cumplen: $x_2 = -x_1^2$
$8x_1^2 + 8x_2 = 8 x_1^2 + 8(-x_1^2) = 0$
por lo que no podemos concluir...
Usemos una función de SciPy scipy.optimize.fmin
```{margin}
Elegimos un punto inicial.
```
End of explanation
"""


print(res_fmin)

print(gf_eval(res_fmin))

Hf_fmin = Hf_eval(res_fmin)

"""
Explanation: ```{margin}
El punto res_fmin satisface $\nabla f_o(x) = 0$ y la relación $x_2 = -x_1^2$.
```
End of explanation
"""


print(np.linalg.eigvals(Hf_fmin))

"""
Explanation: ```{margin}
La Hessiana en res_fmin es semidefinida positiva por lo que no podemos concluir que sea mínimo local de $f_o$.
```
End of explanation
"""


x1_plot,x2_plot = np.meshgrid(np.linspace(-2,2,100), np.linspace(-4,2,100))
z_plot = x1_plot**4 + 2*x1_plot**2*x2_plot + x2_plot**2

point1 = (root1[0], root1[1], fo_numpy(root1))
point2 = (root2[0], root2[1], fo_numpy(root2))
point3 = (root3[0], root3[1], fo_numpy(root3))
#another point:
point4 = (-2, -4, 0)

x1_plot2 = np.linspace(-2,2,100)

# Create the figure
fig = plt.figure()
# Add an axes
ax = fig.gca(projection='3d')
ax.plot_surface(x1_plot, x2_plot, z_plot, alpha=0.2)
ax.scatter(point[0], point[1], point[2], color='green')
ax.scatter(point1[0], point1[1], point1[2], color='green')
ax.scatter(point2[0], point2[1], point2[2], color='green')
ax.scatter(point3[0], point3[1], point3[2], color='green')
ax.scatter(point4[0], point4[1], point4[2], color='green')
ax.plot(x1_plot2, -x1_plot2**2, color="red")
plt.title("$f_o(x) = x_1^4+2x_1^2x_2+x_2^2$")
plt.show()

"""
Explanation: Grafiquemos los puntos que cumplen $x_2=-x_1^2$
End of explanation
"""


x1_plot,x2_plot = np.meshgrid(np.linspace(-100,100,100), np.linspace(-100,100,100))
z_plot = x1_plot**4 + 2*x1_plot**2*x2_plot + x2_plot**2

# Create the figure
fig = plt.figure()
# Add an axes
ax = fig.gca(projection='3d')
ax.plot_surface(x1_plot, x2_plot, z_plot, alpha=0.2)
plt.title("$f_o(x) = x_1^4+2x_1^2x_2+x_2^2$")
plt.show()

np.sum(z_plot < 0)

"""
Explanation: Evaluemos en un rango más amplio la función objetivo $f_o$ y realicemos una gráfica
End of explanation
"""


pointx = np.array([-.5, -1.5])
pointy = np.array([.5, -1.5])
theta = 1/2

point_convex_combination = theta*pointx + (1-theta)*pointy

print(fo_numpy(point_convex_combination))

"""
Explanation: El mínimo valor de $f_o$ es $0$ por lo que tenemos un conjunto de mínimos dados por la curva $x_2=-x_1^2$.
¿Por qué fue un poco complicado determinar el conjunto de mínimos de $f_o$?
$f_o(x)=x_1^4+2x_1^2x_2+x_2^2$ no es una función convexa en su $\text{dom}f_o$, $\mathbb{R}^2$, esto es, no satisface la desigualdad:
$$f_o(\theta x + (1-\theta) y) \leq \theta f_o(x) + (1-\theta)f_o(y)$$
$\forall x,y$ en su dominio y $\forall \theta \in [0,1]$.
End of explanation
"""


print(theta*fo_numpy(pointx) + (1-theta)*fo_numpy(pointy))

"""
Explanation: ```{margin}
Se tiene $f(\theta x + (1-\theta)y) > \theta f(x) + (1-\theta)f(y)$ con la elección $\theta=0.5$, $x=[-0.5, -1.5]^T$, $y=[0.5, -1.5]^T$ lo que indica que $f$ no es convexa sino cóncava para esta elección de puntos $x,y$.
```
End of explanation
"""


P = sympy.Matrix([[5, 4],
                  [4, 5]])

x = sympy.Matrix(sympy.symbols("x1, x2"))

q = sympy.Matrix([-1,1])


r = 3

fo_sympy = (1/2*x.T*P*x + q.T*x)[0] + r

sympy.pprint(fo_sympy.expand())

"""
Explanation: ```{admonition} Observación
:class: tip
Recordar que si $f_o$ es una función convexa, el gradiente de $f_o$ nos ayuda a determinar si un punto es un mínimo local de forma necesaria y suficiente.
```
```{admonition} Ejercicio
:class: tip
Realizar un análisis similar al anterior para la función $f_o(x) = x_2^4+2x_2^2x_1^2+x_1^2$.
```
Ejemplo importante
$$\displaystyle \min_{x \in \mathbb{R}^2} \frac{1}{2}x^TPx+q^Tx+r$$
donde: $P=\left [\begin{array}{cc} 5 & 4\ 4 & 5 \end{array} \right ]$, $q=\left [\begin{array}{c} -1\ 1 \end{array} \right]$, $r=3$.
Haciendo las multiplicaciones de matriz-vector y productos punto
End of explanation
"""


sympy.pprint(P.eigenvals())

"""
Explanation: se reescribe el problema de optimización como:
$$\displaystyle \min_{x \in \mathbb{R}^2} \frac{5}{2}x_1^2 + \frac{5}{2}x_2^2+4x_1x_2 -x_1 + x_2+3$$ 
La función objetivo es una función estrictamente convexa en $\mathbb{R}^2$ (de hecho fuertemente convexa) pues:
```{margin}
Los eigenvalores de $P$ son $1$ y $9$ de multiplicidad simple cada uno.
```
End of explanation
"""


sympy.pprint(derive_by_array(derive_by_array(fo_sympy, (x1,x2)),
                             (x1,x2))
             )

"""
Explanation: y la Hessiana de $f_o$ es:
End of explanation
"""


sympy.pprint(fo_sympy.diff(x))

"""
Explanation: El gradiente de $f_o$ es:
End of explanation
"""


P = np.array([[5,4],[4,5]])
q = np.array([-1,1])
print(np.linalg.solve(P,-q))

"""
Explanation: El mínimo debe satisfacer la ecuación lineal con dos variables y dos ecuaciones
$$\nabla f_o(x) = Px +q =\left [ \begin{array}{cc}
5 & 4\
4 & 5
\end{array}
\right ]
\left [ \begin{array}{c}
x_1\
x_2
\end{array}
\right ]
+ \left [ \begin{array}{c}
-1\
1
\end{array}
\right ]=
\left [ \begin{array}{cc}
5x_1+4x_2-1\
4x_1+5x_2+1
\end{array}
\right ]
=0
$$
```{admonition} Observación
:class: tip
En algunos casos especiales es posible resolver la ecuación no lineal $\nabla f_o(x) = 0$ para $x$ de forma analítica o cerrada. Este es el caso de este ejemplo cuya solución está dada por $x^* = -P^{-1}q$.
```
End of explanation
"""


n = 2 #number of variables
x = cp.Variable(n) #optimization variable
fo_cvxpy = (1/2)*cp.quad_form(x, P) + q.T @ x + r#objective function

opt_objective = cp.Minimize(fo_cvxpy) #optimization objective

prob = cp.Problem(opt_objective) #optimization problem

print(prob.solve())

print("status:", prob.status)
print("optimal value", prob.value)
print("optimal var", x.value)

"""
Explanation: El problema anterior también lo podemos resolver con cvxpy pues es un problema convexo.
```{margin}
Definición de variables y función objetivo: $\frac{1}{2}x^TPx+q^Tx+r$
```
End of explanation
"""


x1, x2, x3, x4 = sympy.symbols("x1, x2, x3, x4")

fo_sympy = (x1-2)**2 + (2-x2)**2 + x3**2 + x4**4

gf = derive_by_array(fo_sympy, (x1, x2, x3, x4))

sympy.pprint(gf)

Hf = derive_by_array(gf, (x1, x2, x3, x4))

sympy.pprint(Hf)

"""
Explanation: ````{admonition} Observaciones
:class: tip

El paquete cvxpy requiere que se especifique el problema de optimización a resolver siguiendo reglas establecidas en Disciplined Convex Programming. En el ejemplo anterior resulta en un error si se ejecutan las líneas siguientes:

```python
x1 = cp.Variable()
x2 = cp.Variable()
fo = 2.5x12 + 4x1x2 - x1 + 2.5x2**2 + x2 + 3
obj = cp.Minimize(fo)
prob = cp.Problem(obj)
prob.solve()
```
La última línea produce
Problem does not follow DCP rules. Specifically: The objective is not DCP.

En la liga de ejemplos hay muchos problemas típicos en optimización convexa y en Atomic Functions ejemplos de funciones atómicas que pueden aplicarse a expresiones de CVXPY.

````
```{admonition} Ejercicio
:class: tip
Utilizando las herramientas (teóricas y prácticas) del ejemplo anterior, resolver el problema de optimización:
$$\min_{x \in \mathbb{R}^2} ||Ax-b||_2^2$$
con $A=\left [ \begin{array}{ccc} 1 & 6 & 2.5\ 1 & 2 & 8 \ 1 & 10 & -1\ 1 & -9 & 3\ 1 & -1 & 2 \end{array} \right ]$, $b=\left [ \begin{array}{c} -1 \ 0 \ 2 \ 3.5 \ -1.7 \end{array} \right ]$.
```
Métodos de descenso
Los ejemplos anteriores mostraron la importancia de la información de primer y segundo orden de la función objetivo $f_o$ y las ecuaciones no lineales para resolver el problema de optimización. El primer ejemplo utilizó las funciones scipy.optimize.fsolve y scipy.optimize.fmin para este propósito. Tales funciones utilizan métodos iterativos para resolver ecuaciones no lineales y aproximar un mínimo local respectivamente en un esquema que para el caso de minimización satisface:
$$f_o(x^{(0)}) > f_o(x^{(1)}) > f_o(x^{(2)}) > \cdots > f_o(x^{(k)}) > \cdots$$
con $x^{(0)}$ punto inicial.
Los métodos que utilizan esquemas iterativos para calcular una secuencia de minimización de puntos $x^{(0)}, x^{(1)}, \dots \in \text{dom}f_o$ con la característica $f_o(x^{(k)}) \rightarrow p^$ si $k \rightarrow \infty$ se conocen con el nombre de métodos de descenso*.
```{admonition} Definición
Un método de descenso es aquel que genera la secuencia de minimización $x^{(0)}, x^{(1)}, \dots \in \text{dom}f_o$ la cual cumple con la desigualdad: $f_o(x^{(k+1)}) < f_o(x^{(k)})$ excepto para $x^{(k)}$ óptimo y $f_o(x^{(k)}) \rightarrow p^*$ si $k \rightarrow \infty$.
```
La secuencia de minimización en un método de descenso se genera con la expresión:
$$x^{(k+1)} = x^{(k)} + \Delta x$$
donde: $\Delta x$ es dirección de búsqueda que es de descenso.
```{admonition} Observaciones
:class: tip


La definición de una dirección de descenso anterior aplica para funciones en general, no es necesario que $f_o$ sea convexa.


En la práctica la secuencia de minimización involucra calcular un tamaño o longitud de paso que ayuda a la convergencia del método. Se tienen diferentes metodologías para su cálculo y más adelante se describe una de ellas.


```
```{margin}
Recuérdese que si $f_o$ es fuertemente convexa en el conjunto $\mathcal{S}$ entonces $\nabla^2 f_o (x) \in \mathbb{S}^n_{++}$ y $\text{cond}(\nabla ^2 f_o(x))$ está acotado por arriba por una constante para $x \in \mathcal{S}$.
```
```{admonition} Comentario
Asumiendo que la función $f_o$ es convexa, típicamente se asume lo siguiente para tener métodos iterativos confiables y exactos:


Los puntos iniciales $x^{(0)}$ están en $\text{dom}f_o$.


Que el conjunto $f_o(x^{(0)})$-subnivel sea cerrado pues así se garantiza que la secuencia de minimización está en el conjunto $f_o(x^{(0)})$-subnivel para todas las iteraciones.


$f_o$ fuertemente convexa en el conjunto $f_o(x^{(0)})$-subnivel para tener propiedades dadas en los {ref}resultados que son posibles probar para funciones fuertemente convexas &lt;RESFFUERTCON&gt;. 


```
En lo siguiente se asume que $f_o$ cumple $f_o \in \mathcal{C}^2(\text{dom}f_o)$ y es convexa en un conjunto convexo y cerrado que contiene a $x^*$. Ver {ref}conjunto abierto, cerrado, cerradura e interior &lt;CACCI&gt; para definición de conjunto cerrado.
Ejemplo de función objetivo convexa
Encontrar el mínimo del siguiente problema con un método iterativo.
$$\displaystyle \min_{x \in \mathbb{R}^4} (x_1-2)^2+(2-x_2)^2+x_3^2+x_4^4$$
End of explanation
"""


x_0 = np.array([5,5,1,0.1])

"""
Explanation: Como $f_o$ es una función convexa (de hecho estrictamente convexa) en su dominio $\mathbb{R}^4$, se tiene que su óptimo se obtiene igualando y resolviendo la ecuación no lineal $\nabla f_o(x) = 0$ :
$$\nabla f_o(x) = 
\left[ \begin{array}{c}
2x_1-4 \
2x_2-4\
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 con un método iterativo?
Condición para que un paso o dirección de búsqueda sea de descenso
La idea de los métodos de optimización es calcular direcciones $\Delta x$ de búsqueda que sean de descenso, esto es, que al movernos de un punto a otro en tal dirección, el valor de $f_o$ decrece. Existen muchas direcciones de descenso (de hecho infinitas) una que se muestra en el dibujo siguiente es la dirección de descenso de Newton $\Delta x_{nt}$: 
<img src="https://dl.dropboxusercontent.com/s/25bmebx645howjw/direccion_de_descenso_de_Newton_1d.png?dl=0" heigth="600" width="600">
En el dibujo $f = f_o$ y $\hat{f}$ es un modelo cuadrático. Del punto $(x,f(x))$ nos debemos mover al punto $(x+\Delta x_{nt}, f(x + \Delta x_{nt}))$ para llegar al óptimo. En tal dirección $f$ decrece: $f(x+\Delta x_{nt}) < f(x)$ y obsérvese que $\Delta x_{nt}$ es mínimo de $\hat{f}$.
```{margin}
Ver {ref}teorema de Taylor para una función de varias variables &lt;TEOTAYLORNVARIABLES&gt;
```
```{admonition} Comentario
El modelo cuadrático del dibujo anterior está dado por la aproximación de segundo orden a la función $f_o$ por el teorema de Taylor con centro en $x$:
$$m(x + v) = \hat{f}_o(x + v) = f_o(x) + \nabla f_o(x)^T v + \frac{1}{2} v^T \nabla^2f_o(x)v$$
con único mínimo si $\nabla ^2 f_o(x) \in \mathbb{S}^n_{++}$ dado por $v^* = \Delta x_{nt}$ y $\Delta x_{nt}$ dirección de Newton cuya expresión está más adelante. En cada iteración se construye un modelo cuadrático:
$$m(x^{(k)} + v) = \hat{f}_o(x^{(k)} + v) = f_o(x^{(k)}) + \nabla f_o(x^{(k)})^T v + \frac{1}{2} v^T \nabla^2f_o(x^{(k)})v$$
```
Geométricamente las direcciones de descenso forman un ángulo agudo con $-\nabla f_o(x)$:
<img src="https://dl.dropboxusercontent.com/s/eednhn6lj1rag1j/zone_for_descent_directions.png?dl=0" heigth="350" width="350">
En el dibujo $f = f_o$.
```{admonition} Observación
:class: tip
Aunque se tienen una cantidad infinita de direcciones de descenso, las direcciones de descenso que típicamente son elegidas no son cercanas a ser ortogonales con el gradiente de $f_o$.
```
Tenemos una condición para garantizar que una dirección sea de descenso:
```{admonition} Definición
Si el paso o dirección de búsqueda satisface: $\nabla f_o(x)^T\Delta x < 0$ se le nombra dirección de descenso. 
```
```{margin}
Recuérdese que el teorema de Taylor nos ayuda a aproximar a una función de forma local.
```
```{admonition} Comentarios


Recuérdese que $\nabla f_o(x)^T \Delta x$ es una derivada direccional de $f_o$ en $x$ en la dirección $\Delta x$, ver {ref}ejemplo función restringida a una línea &lt;EJRestriccionALinea&gt;.


La definición anterior se justifica pues recuérdese que por la aproximación del teorema de Taylor a primer orden se tiene:


$$f_o(x + \Delta x) \approx f_o(x) + \nabla f_o(x) ^T \Delta x$$
y si $\Delta x$ es dirección de descenso entonces: $f_o(x) + \nabla f_o(x) ^T \Delta x < f_o(x)$.
```
```{admonition} Observación
:class: tip
Obsérvese que si $x^$ es mínimo local entonces $\nabla f_o(x^) = 0$ (condición necesaria de primer orden) por lo que no existen direcciones de descenso.
```
Ejemplos de direcciones de descenso
```{sidebar} La dirección del gradiente...
Una forma de obtener la dirección de Newton es encontrando el mínimo del modelo cuadrático referido anteriormente asumiendo $\nabla ^2 f_o(x) \in \mathbb{S}^n_{++}$. La dirección del gradiente se obtiene al resolver el problema de optimización con restricciones siguiente:
$$\min_{v \in \mathbb{R}^n} \nabla f_o(x)^T v$$
$$\text{sujeto a:} ||v|| = 1$$
para la norma $2$. Se utiliza una restricción del tipo normalización pues la función objetivo involucra un producto punto (derivada direccional) que es dependiente o proporcional a la longitud de $v$.
```


$\Delta x = - \nabla f_o \left (x^{(k)} \right )$ que da lugar al método de descenso en gradiente para $x^{(k)}$ no óptimo.


$\Delta x = - \nabla^2 f_o \left (x^{(k)} \right )^{-1} \nabla f_o\left(x^{(k)} \right)$ que da lugar al método de descenso por Newton con $\nabla^2 f_o \left (x^{(k)} \right ) \in \mathbb{S}^n_{++}$ y $x^{(k)}$ no óptimo.


$\Delta x = - H_k ^{-1} \nabla f_o\left(x^{(k)}\right)$  con $H_k$ aproximación a la Hessiana de $f_o$ con $\nabla^2 f_o \left (x^{(k)} \right ) \in \mathbb{S}^n_{++}$ y $x^{(k)}$ no óptimo.


```{admonition} Observaciones
:class: tip


La definición de una dirección de descenso aplica para funciones en general, no es necesario que $f_o$ sea convexa.


Para funciones en general, la dirección de Newton es de descenso si la Hessiana es definida positiva y análogamente para las direcciones en las que se utilicen aproximaciones a la Hessiana. Esto asegura que el modelo cuadrático tenga un único mínimo y que $f_o$ decrezca su valor en tal dirección.


Comúnmente los métodos que utilizan aproximaciones a la Hessiana se conocen con el nombre de métodos Cuasi-Newton, ver Quasi-Newton_method.
```


Continuando con el ejemplo anterior
Encontrar el mínimo del siguiente problema con un método iterativo.
$$\displaystyle \min_{x \in \mathbb{R}^4} (x_1-2)^2+(2-x_2)^2+x_3^2+x_4^4$$
Opción descenso en gradiente: usando la dirección del gradiente de $f_o$ se tiene: 
$$x^{(k+1)} = x^{(k)} - \nabla f_o(x^{(k)})$$
Tomando $x^{(0)} = \left [ \begin{array}{c} 5 \ 5 \1 \ 0.1 \ \end{array} \right ]$ como punto inicial:
End of explanation
"""


f_o_np = lambda x: (x[0]-2)**2 + (2-x[1])**2 + x[2]**2 + x[3]**2

"""
Explanation: ```{margin}
Definición de $f_o$
```
End of explanation
"""


print(f_o_np(x_0))

gf_eval = lambda x: np.array([partial_derivative.subs({"x1": x[0], 
                                                       "x2": x[1],
                                                       "x3": x[2],
                                                       "x4": x[3]}) for partial_derivative in gf],
                              dtype=float)

"""
Explanation: ```{margin}
Evaluando $f_o$ en $x^{(0)}$.
```
End of explanation
"""


print(gf_eval(x_0))

"""
Explanation: ```{margin}
Evaluando el gradiente en $x^{(0)}$: $\nabla f_o(x^{(0)})$.
```
End of explanation
"""


print(gf_eval(x_0).dot(-gf_eval(x_0)))

"""
Explanation: ```{margin}
Verificando que es dirección de descenso: $\nabla f_o \left (x^{(0)} \right )^T \Delta x < 0$ con $\Delta x = -\nabla f_o \left (x^{(0)} \right )$.
```
End of explanation
"""


x_1 = x_0 - gf_eval(x_0)

print(x_1)

"""
Explanation: Primera iteración
```{margin}
Esquema iterativo: $x_1 = x_0 + \Delta x = x_0 - \nabla f_o(x^{(0)})$.
```
End of explanation
"""


print(f_o_np(x_0))

print(f_o_np(x_1))

"""
Explanation: ```{margin}
Evaluando $f_o$ en $x^{(1)}$ se observa que $f_o$ decreció muy poco, de hecho $f_o(x^{(0)}) \approx f_o(x^{(1)})$.
```
End of explanation
"""


print(f_o_np(x_0) + gf_eval(x_0).dot(-gf_eval(x_0)))

"""
Explanation: Obsérvese que la aproximación a primer orden por Taylor no es correcta: $f_o(x_0 + \Delta x) \neq f_o(x_0) + \nabla f_o(x_0)^T \Delta x$.
End of explanation
"""


print(gf_eval(x_1))

"""
Explanation: Segunda iteración
```{margin}
Evaluando el gradiente en $x^{(1)}$: $\nabla f_o(x^{(1)})$.
```
End of explanation
"""


print(gf_eval(x_1).dot(-gf_eval(x_1)))

"""
Explanation: ```{margin}
Verificando que es dirección de descenso: $\nabla f_o \left (x^{(1)} \right )^T \Delta x < 0$ con $\Delta x = -\nabla f_o \left (x^{(1)} \right )$.
```
End of explanation
"""


x_2 = x_1 - gf_eval(x_1)

print(x_2)

"""
Explanation: ```{margin}
Esquema iterativo: $x_2 = x_1 + \Delta x = x_1 - \nabla f_o(x^{(1)})$.
```
End of explanation
"""


print(f_o_np(x_1))

print(f_o_np(x_2))

"""
Explanation: ```{margin}
Evaluando $f_o$ en $x^{(1)}$ se observa que $f_o$ decreció muy poco, de hecho $f_o(x^{(2)}) \approx f_o(x^{(1)})$.
```
End of explanation
"""


print(f_o_np(x_1) + gf_eval(x_1).dot(-gf_eval(x_1)))

"""
Explanation: Obsérvese que la aproximación a primer orden por Taylor no es correcta: $f_o(x_1 + \Delta x) \neq f_o(x_1) + \nabla f_o(x_1)^T \Delta x$.
End of explanation
"""


x_3 = x_2 - gf_eval(x_2)

"""
Explanation: Tercera iteración
```{margin}
Esquema iterativo: $x_3 = x_2 + \Delta x = x_2 - \nabla f_o(x^{(2)})$.
```
End of explanation
"""


print(x_1)

print(x_3)

"""
Explanation: ```{margin}
Obsérvese que $x_3 \approx x_1$.
```
End of explanation
"""


x_4 = x_3 - gf_eval(x_3)

"""
Explanation: Cuarta iteración
```{margin}
Esquema iterativo: $x_4 = x_3 + \Delta x = x_3 - \nabla f_o(x^{(3)})$.
```
End of explanation
"""


print(x_2)

print(x_4)

"""
Explanation: ```{margin}
Obsérvese que $x_4 \approx x_2$.
```
End of explanation
"""


x_0 = np.array([5,5,1,0.1])

Hf_eval = lambda x: np.array([second_partial_derivative.subs({"x1": x[0], 
                                                              "x2": x[1],
                                                              "x3": x[2],
                                                              "x4": x[3]}) for second_partial_derivative in Hf],
                             dtype=float)

"""
Explanation: y así nos quedaremos ciclando por muchas iteraciones...
```{admonition} Comentario
El método de descenso en gradiente para el ejemplo anterior no convergerá al óptimo $x^* = \left [ \begin{array}{c} 2 \ 2 \0 \ 0 \ \end{array} \right ]$
```
```{admonition} Ejercicio
:class: tip
¿Converge el método de descenso en gradiente al óptimo para un punto inicial $x^{(0)} = \left [ \begin{array}{c} 2.5 \ 2.5 \0.001 \ 0.001 \ \end{array} \right ]$?
```
Opción descenso por dirección de Newton: usando la dirección de descenso de Newton de $f_o$ se tiene: 
$$x^{(k+1)} = x^{(k)} - \nabla^2 f_o \left (x^{(k)} \right )^{-1} \nabla f_o\left(x^{(k)} \right)$$
Con 
$$
\nabla ^2 f_o(x) = 
\left [ 
\begin{array}{cccc}
2 & 0 & 0 & 0 \
0 & 2 & 0 & 0 \
0 & 0 & 2 & 0 \
0 & 0 & 0 & 12x_4^2
\end{array}
\right ]
$$
Tomando $x^{(0)} = \left [ \begin{array}{c} 5 \ 5 \1 \ 0.1 \ \end{array} \right ]$ como punto inicial y no calculando la inversa de la Hessiana pues en su lugar resolvemos el sistema de ecuaciones lineales $\nabla ^2 f_o \left ( x^{(k)} \right ) \Delta x = - \nabla f_o \left(x^{(k)}\right )$ resulta
End of explanation
"""


Hf_sympy_eval = sympy.Matrix(Hf_eval(x_0))
sympy.pprint(Hf_sympy_eval)

"""
Explanation: ```{margin}
Evaluando la Hessiana en $x^{(0)}$: $\nabla^2f_o(x^{(0)})$.
```
End of explanation
"""


gf_sympy_eval = gf_eval(x_0)
sympy.pprint(gf_sympy_eval)

"""
Explanation: ```{margin}
Evaluando el gradiente en $x^{(0)}$: $\nabla f_o(x^{(0)})$.
```
End of explanation
"""


Hf_np_eval = np.array(Hf_sympy_eval, dtype=float)

_, n = Hf_np_eval.shape

"""
Explanation: El sistema de ecuaciones lineales a resolver es:
$$
\left [
\begin{array}{cccc}
2 & 0 & 0 & 0 \
0 & 2 & 0 & 0 \
0 & 0 & 2 & 0 \
0 & 0 & 0 & 0.12
\end{array}
\right ]
\Delta x = 
-\left [
\begin{array}{c}
6 \
6 \
2 \
0.004
\end{array}
\right ]
$$
Resolviendo con NumPy el sistema de ecuaciones lineales:
```{margin} 
Aquí convertimos de un objeto de SymPy a un array de NumPy pues si bien podríamos resolver el sistema con SymPy es menos costoso utilizar arrays de NumPy.
```
End of explanation
"""


gf_np_eval = np.array(gf_sympy_eval, dtype = float)

gf_np_eval

"""
Explanation: ```{margin} 
Aquí convertimos de un objeto de SymPy a un array de NumPy pues si bien podríamos resolver el sistema con SymPy es menos costoso utilizar arrays de NumPy.
```
End of explanation
"""


dir_Newton = np.linalg.solve(Hf_np_eval, -gf_np_eval)

print(dir_Newton)

"""
Explanation: ```{margin}
Resolvemos $\nabla ^2 f_o \left ( x^{(0)} \right ) \Delta x = - \nabla f_o \left(x^{(0)}\right )$ para obtener $\Delta x$ dirección de Newton.
```
End of explanation
"""


print(gf_np_eval.dot(dir_Newton))

"""
Explanation: ```{margin}
Verificando que es dirección de descenso: $\nabla f_o \left (x^{(0)} \right )^T \Delta x < 0$ con $\Delta x = -\nabla ^2f_o(x^{(0)})^{-1} \nabla f_o(x^{(0)})$.
```
End of explanation
"""


x_1 = x_0 + dir_Newton

print(x_1)

"""
Explanation: Primera iteración
```{margin}
Esquema iterativo: $x_1 = x_0 + \Delta x = x_0 - \nabla ^2f_o(x^{(0)})^{-1} \nabla f_o(x^{(0)})$
```
End of explanation
"""


print(np.linalg.cond(Hf_np_eval))

"""
Explanation: Recuérdese que siempre es útil monitorear el número de condición de la matriz del sistema de ecuaciones lineales que en este caso es la Hessiana de $f_o$ en $x^{(0)}$ para confiabilidad de nuestros cálculos al resolver el sistema de ecuaciones lineales asociado, ver {ref}Número de condición de una matriz &lt;NCM&gt;:
End of explanation
"""


print(f_o_np(x_0))

print(f_o_np(x_1))

"""
Explanation: ```{margin}
Evaluando $f_o$ en $x^{(1)}$ se observa que $f_o$ sí decrece $f_o(x^{(1)}) < f_o(x^{(0)})$.
```
End of explanation
"""


print(f_o_np(x_0) + gf_np_eval.dot(dir_Newton) + 1/2*dir_Newton.dot(Hf_np_eval@dir_Newton))

"""
Explanation: Obsérvese que la aproximación a segundo orden por Taylor es cercana: $f_o(x_0 + \Delta x) \approx f_o(x_0) + \nabla f_o(x_0)^T \Delta x + \frac{1}{2}\Delta x ^T\nabla^2 f_o(x_0) \Delta x$.
End of explanation
"""


Hf_sympy_eval = sympy.Matrix(Hf_eval(x_1))
gf_sympy_eval = gf_eval(x_1)
Hf_np_eval = np.array(Hf_sympy_eval, dtype=float)
gf_np_eval = np.array(gf_sympy_eval, dtype = float)
dir_Newton = np.linalg.solve(Hf_np_eval, -gf_np_eval)
x_2 = x_1 + dir_Newton

print(x_2)

"""
Explanation: Segunda iteración
```{margin}
Evaluando la Hessiana en $x^{(1)}$: $\nabla^2f_o(x^{(1)})$, evaluando el gradiente en $x^{(1)}$: $\nabla f_o(x^{(1)})$ y resolviendo el sistema de ecuaciones $\nabla ^2 f_o \left ( x^{(1)} \right ) \Delta x = - \nabla f_o \left(x^{(1)}\right )$.
```
```{margin}
Esquema iterativo: $x_2 =  x_1 + \Delta x = x_1 - \nabla ^2f_o(x^{(1)})^{-1} \nabla f_o(x^{(1)})$
```
End of explanation
"""


print(gf_np_eval.dot(dir_Newton))

"""
Explanation: ```{margin}
Verificando que es dirección de descenso: $\nabla f_o \left (x^{(1)} \right )^T \Delta x < 0$ con $\Delta x = -\nabla ^2f_o(x^{(1)})^{-1} \nabla f_o(x^{(1)})$.
```
End of explanation
"""


print(f_o_np(x_1))

print(f_o_np(x_2))

"""
Explanation: ```{margin}
Evaluando $f_o$ en $x^{(2)}$ se observa que $f_o$ sí decrece $f_o(x^{(2)}) < f_o(x^{(1)})$.
```
End of explanation
"""


print(f_o_np(x_1) + gf_np_eval.dot(dir_Newton) + 1/2*dir_Newton.dot(Hf_np_eval@dir_Newton))

"""
Explanation: Obsérvese que la aproximación a segundo orden por Taylor es cercana: $f_o(x_1 + \Delta x) \approx f_o(x_1) + \nabla f_o(x_1)^T \Delta x + \frac{1}{2}\Delta x ^T\nabla^2 f_o(x_1) \Delta x$.
End of explanation
"""


print(np.linalg.cond(Hf_np_eval))

"""
Explanation: ```{margin}
Número de condición de la Hessiana de $f_o$ en $x^{(1)}$.
```
End of explanation
"""


Hf_sympy_eval = sympy.Matrix(Hf_eval(x_2))
gf_sympy_eval = gf_eval(x_2)
Hf_np_eval = np.array(Hf_sympy_eval, dtype=float)
gf_np_eval = np.array(gf_sympy_eval, dtype = float)
dir_Newton = np.linalg.solve(Hf_np_eval, -gf_np_eval)
x_3 = x_2 + dir_Newton

print(x_3)

"""
Explanation: Tercera iteración
```{margin}
Evaluando la Hessiana en $x^{(2)}$: $\nabla^2f_o(x^{(2)})$, evaluando el gradiente en $x^{(2)}$: $\nabla f_o(x^{(2)})$  y resolviendo el sistema de ecuaciones $\nabla ^2 f_o \left ( x^{(2)} \right ) \Delta x = - \nabla f_o \left(x^{(2)}\right )$.
```
```{margin}
Esquema iterativo: $x_3 = x_2 + \Delta x = x_2 - \nabla ^2f_o(x^{(2)})^{-1} \nabla f_o(x^{(2)})$
```
End of explanation
"""


print(gf_np_eval.dot(dir_Newton))

"""
Explanation: ```{margin}
Verificando que es dirección de descenso: $\nabla f_o \left (x^{(2)} \right )^T \Delta x < 0$ con $\Delta x = -\nabla ^2f_o(x^{(2)})^{-1} \nabla f_o(x^{(2)})$.
```
End of explanation
"""


print(f_o_np(x_2))

print(f_o_np(x_3))

"""
Explanation: ```{margin}
Evaluando $f_o$ en $x^{(3)}$ se observa que $f_o$ sí decrece $f_o(x^{(3)}) < f_o(x^{(2)})$.
```
End of explanation
"""


print(f_o_np(x_2) + gf_np_eval.dot(dir_Newton) + 1/2*dir_Newton.dot(Hf_np_eval@dir_Newton))

"""
Explanation: Obsérvese que la aproximación a segundo orden por Taylor es cercana: $f_o(x_2 + \Delta x) \approx f_o(x_2) + \nabla f_o(x_2)^T \Delta x + \frac{1}{2}\Delta x ^T\nabla^2 f_o(x_2) \Delta x$.
End of explanation
"""


print(np.linalg.cond(Hf_np_eval))

"""
Explanation: ```{margin}
Número de condición de la Hessiana de $f_o$ en $x^{(2)}$.
```
End of explanation
"""


Hf_sympy_eval = sympy.Matrix(Hf_eval(x_3))
gf_sympy_eval = gf_eval(x_3)
Hf_np_eval = np.array(Hf_sympy_eval, dtype=float)
gf_np_eval = np.array(gf_sympy_eval, dtype = float)
dir_Newton = np.linalg.solve(Hf_np_eval, -gf_np_eval)
x_4 = x_3 + dir_Newton

print(x_4)

"""
Explanation: Cuarta iteración
```{margin}
Evaluando la Hessiana en $x^{(3)}$: $\nabla^2f_o(x^{(3)})$, evaluando el gradiente en $x^{(3)}$: $\nabla f_o(x^{(3)})$  y resolviendo el sistema de ecuaciones $\nabla ^2 f_o \left ( x^{(3)} \right ) \Delta x = - \nabla f_o \left(x^{(3)}\right )$.
```
```{margin}
Esquema iterativo: $x_4 = x_3 + \Delta x = x_3 - \nabla ^2f_o(x^{(3)})^{-1} \nabla f_o(x^{(3)})$
```
End of explanation
"""


print(gf_np_eval.dot(dir_Newton))

"""
Explanation: ```{margin}
Verificando que es dirección de descenso: $\nabla f_o \left (x^{(3)} \right )^T \Delta x < 0$ con $\Delta x = -\nabla ^2f_o(x^{(3)})^{-1} \nabla f_o(x^{(3)})$.
```
```{margin}
Obsérvese que el gradiente y la dirección de Newton son cercanos a ser ortogonales.
```
End of explanation
"""


print(f_o_np(x_3))

print(f_o_np(x_4))

"""
Explanation: ```{margin}
Evaluando $f_o$ en $x^{(4)}$ se observa que $f_o$ sí decrece $f_o(x^{(4)}) < f_o(x^{(3)})$.
```
End of explanation
"""


print(f_o_np(x_3) + gf_np_eval.dot(dir_Newton) + 1/2*dir_Newton.dot(Hf_np_eval@dir_Newton))

"""
Explanation: Obsérvese que la aproximación a segundo orden por Taylor es cercana: $f_o(x_3 + \Delta x) \approx f_o(x_3) + \nabla f_o(x_3)^T \Delta x + \frac{1}{2}\Delta x ^T\nabla^2 f_o(x_3) \Delta x$.
End of explanation
"""


print(np.linalg.cond(Hf_np_eval))

print(x_4 == approx(np.array([2,2,0,0.0]), abs=1e-1, rel=1e-1))

"""
Explanation: ```{margin}
Obsérvese cómo va aumentando el número de condición de la Hessiana conforme nos aproximamos a la solución, en este paso se ha calculado la Hessiana de $f_o$ en $x^{(3)}$.
```
End of explanation
"""


x_0 = np.array([5, 5, 1, 0])

"""
Explanation: ```{admonition} Comentario
El método por dirección de Newton sí convergerá al óptimo $x^* = \left [ \begin{array}{c} 2 \ 2 \0 \ 0 \ \end{array} \right ]$ pero la convergencia será lenta.
```
Si hubiéramos elegido como punto inicial $x^{(0)} = \left [ \begin{array}{c} 5 \ 5 \1 \ 0 \ \end{array} \right ]$
End of explanation
"""


Hf_sympy_eval = sympy.Matrix(Hf_eval(x_0))
sympy.pprint(Hf_sympy_eval)

"""
Explanation: ```{margin}
Evaluando la Hessiana en $x^{(0)}$: $\nabla^2f_o(x^{(0)})$.
```
End of explanation
"""


gf_sympy_eval = sympy.Matrix(gf_eval(x_0))
sympy.pprint(gf_sympy_eval)

"""
Explanation: ```{margin}
Evaluando el gradiente en $x^{(0)}$: $\nabla f_o(x^{(0)})$.
```
End of explanation
"""


x = sympy.Matrix([x1, x2, x3, x4])

"""
Explanation: El sistema de ecuaciones lineales a resolver es:
$$
\left [
\begin{array}{cccc}
2 & 0 & 0 & 0 \
0 & 2 & 0 & 0 \
0 & 0 & 2 & 0 \
0 & 0 & 0 & 0
\end{array}
\right ]
\Delta x = 
-\left [
\begin{array}{c}
6 \
6 \
2 \
0
\end{array}
\right ]
$$
Obsérvese que la última ecuación es redundante por lo que una solución al sistema anterior es considerar $x_4=0$ y resolver:
$$
\left [
\begin{array}{ccc}
2 & 0 & 0  \
0 & 2 & 0  \
0 & 0 & 2  \
\end{array}
\right ]
\Delta x = 
-\left [
\begin{array}{c}
6 \
6 \
2 \
\end{array}
\right ]
$$
End of explanation
"""


sympy.pprint(sympy.solve(Hf_sympy_eval*x - (-gf_sympy_eval), x))

"""
Explanation: ```{margin}
SymPy nos permite obtener soluciones a sistemas de ecuaciones lineales que tienen un renglón y columna de ceros.
```
End of explanation
"""


x_1 = x_0 + np.array([-3, -3, -1, 0])

"""
Explanation: ```{margin}
Esquema iterativo: $x_1 = x_0 + \Delta x = x_0 - \nabla ^2f_o(x^{(0)})^{-1} \nabla f_o(x^{(0)})$
```
End of explanation
"""


print(x_1)

"""
Explanation: ```{margin}
$x_1$ es el óptimo del problema.
```
End of explanation
"""


sympy.pprint(Hf_sympy_eval.condition_number())

print(x_1 == approx(np.array([2,2,0,0.0]), abs=1e-4, rel=1e-4))

"""
Explanation: ```{margin}
El número de condición es $\infty$, con el print de SymPy se ve así.
```
End of explanation
"""


x1 = cp.Variable()
x2 = cp.Variable()
x3 = cp.Variable()
x4 = cp.Variable()

fo_cvxpy = (x1 -2)**2 + (2-x2)**2 + x3**2 + x4**2



obj = cp.Minimize(fo_cvxpy)

prob = cp.Problem(obj)

print(prob.solve())

print("status:", prob.status)
print("optimal value", prob.value)
print("optimal var", x1.value, x2.value, x3.value, x4.value)

"""
Explanation: ```{admonition} Comentarios
De acuerdo al ejemplo anterior:


Utilizar información de primer o segundo orden nos ayuda a encontrar óptimo(s) de funciones.


El método de descenso en gradiente no converge para el punto inicial elegido. En el caso del método de Newton sí hay convergencia pero es lenta si el punto inicial tiene en la última entrada un número cercano a $0$.


Una de las ventajas que tiene utilizar la dirección del gradiente vs la dirección de Newton es que el gradiente involucra menos almacenamiento en memoria que el almacenamiento de la Hessiana: $\mathcal{O}(n)$ vs $\mathcal{O}(n^2)$. 


```
Resolviendo el problema con CVXPY
End of explanation
"""


x_0 = np.array([5,5,1,0.1])

t_0=0.5

x_1 = x_0 - t_0*gf_eval(x_0)

print(x_1)

print(x_1 == approx(np.array([2,2,0,0.0]), abs=1e-1, rel=1e-1))

"""
Explanation: Tamaño o longitud de paso
En el ejemplo anterior en el que se aproximó al mínimo del siguiente problema con el método de Newton
$$\displaystyle \min_{x \in \mathbb{R}^4} (x_1-2)^2+(2-x_2)^2+x_3^2+x_4^4$$
se concluyó que tal método converge de forma lenta y el método de descenso en gradiente no converge para el punto inicial elegido en ambos métodos. La pequeña reducción que se obtenía en $f_o$ por cada iteración fue la razón de tal situación en el caso del descenso en gradiente. Una metodología que resuelve la no convergencia del método de descenso en gradiente utiliza el siguiente esquema iterativo:
$$x^{(k+1)} = x^{(k)} + t^{(k)}\Delta x$$
con $t^{(k)}>0$.
```{admonition} Comentario
Para el caso del método de Newton, la convergencia cuadrática depende además de elegir tamaños de pasos adecuados que:


las iteraciones vayan aproximándose a $x^*$,


la función objetivo sea dos veces diferenciable y su Hessiana sea definida positiva en un conjunto abierto que contenga a $x^$ y sea Lipschitz* continua en tal conjunto, ver Lipschitz_continuity , que ayuda a acotar la diferencia entre $f_o$ y el modelo cuadrático $m$.


```
Continuando con el ejemplo anterior
Usando el método de descenso en gradiente de $f_o$ y el esquema iterativo: 
$$x^{(k+1)} = x^{(k)} - t^{(k)} \nabla f_o(x^{(k)})$$
con $t^{(0)} = 0.5$ con punto inicial $x^{(0)} = \left [ \begin{array}{c} 5 \ 5 \1 \ 0.1 \ \end{array} \right ]$ se tiene:
End of explanation
"""


f_o = lambda x:x**2
n = 10

def minimization_sequence():
    for k in range(n):
        yield (-1)**k*(1+2**(-k))

t = np.linspace(-2.3, 2.3, 100)

plt.plot(t, f_o(t))
[plt.scatter(s, f_o(s)) for s in minimization_sequence()]
plt.title("Secuencia de minimización que no converge al mínimo de $f_o$")
plt.annotate('$f(x^{(0)})$',(2, f_o(2)),fontsize=12)
plt.annotate('$f(x^{(1)})$',(-1.5, f_o(-1.5)),fontsize=12)
plt.annotate('$f(x^{(2)})$',(1.25, f_o(1.25)),fontsize=12)
plt.show()

"""
Explanation: por lo que llegamos a aproximar al óptimo en una iteración.
```{admonition} Observación
:class: tip
Para problemas small scale el método de Newton o Cuasi-Newton son muy superiores al método de descenso en gradiente. Sin embargo para problemas large scale el método de descenso en gradiente es utilizado ampliamente en aplicaciones de machine learning.
```
```{admonition} Definición
Al escalar $t^{(k)}$ se le nombra tamaño o longitud de paso y siempre es positivo salvo en el caso en que $x^{(k)}$ sea óptimo.
```
El valor $t^{(k)}$ se calcula con metodologías como búsqueda de línea o regiones de confianza, ver line search, trust region y en esta nota se revisa la búsqueda de línea con backtracking.
¿Por qué funciona lo anterior?
La condición para que una dirección produzca descenso en $f_o$, $f_o(x^{(k+1)}) < f_o(x^{(k)})$, no es la única para que la secuencia de minimización converja a un mínimo de $f_o$. El siguiente ejemplo muestra el comportamiento del ejemplo trabajado antes en el plano.
Considérese $f_o(x) = x^2$, $x^{(0)} = 2$ y los pasos dados por la secuencia: $(-1)^k(1+2^{-k})$. Entonces:
End of explanation
"""


