# Optimization

In general, optimization is the process of finding and selecting the optimal element
from a set of feasible candidates. In mathematical optimization, this problem is usually formulated as
determining the extreme value of a function of a given domain. An extreme value, or an optimal value, can
refer to either the minimum or maximum of the function, depending on the application and the specific
problem.

Optimization is closely related to equation solving because at an optimal value of a function, its
derivative, or gradient in the multivariate case, is zero.

a method to solve optimization problems is to solve for the zeros of the derivative or the gradient and test
the resulting candidates for optimality. This approach is not always feasible though, and often it is required
to take other numerical approaches

Here we restrict our attention to mathematical optimization of real-valued functions, with one or more
dependent variables.

## Importing Modules


```python
%matplotlib inline
from scipy import optimize
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

import matplotlib.pyplot as plt
import numpy as np
from sympy import *
init_printing()
```

A general optimization problem of the type considered here can be formulated as a minimization problem,

$$\min_x f(x )$$, 

- subject to sets of $m$ equality constraints $g(x) = 0$ and 
- $p$ inequality constraints $h(x) \leq 0$. 

Here $f(x)$ is a <span class="mark">real-valued function</span> of $x$, which can be a scalar or a vector $x = (x_0,x_1,\ldots,x_n)^T$ while $g(x)$ and $h(x)$ can be <span class="mark">vector valued functions</span>

$f :R^n \rightarrow R, g:R^n \rightarrow R^m $ and $h :R^n \rightarrow R^p$. 

Note that maximizing $f(x)$ is equivalent to minimizing
$–f (x)$, so without loss of generality it is sufficient to consider only minimization problems.

## Optimization problems are classified depending on the properties of the functions $f (x), g(x),$ and $h(x).$


- the problem is univariate or one dimensional if $x$ is a scalar, $x \in R$, 

- multivariate or multidimensional if $x$ is a vector, $x \in R$. For high-dimensional objective functions, with larger $n$, the
optimization problem is harder and more computationally demanding to solve. 

- If the objective function and the constraints all are linear, the problem is a linear optimization problem, or linear programming problem.

- If either the objective function or the constraints are nonlinear, it is a nonlinear optimization problem, or
nonlinear programming problem. 

With respect to constraints, important subclasses of optimization are unconstrained problems, and those with linear and nonlinear constraints. 

Finally, handling equality and inequality constraints require different approaches.

However, an important subclass of nonlinear problems that can be solved efficiently is convex problems, which are directly related to the absence of strictly local minima and the existence of a
unique global minimum. By definition, a function is convex on an interval $[a, b]$ if the values of the function
on this interval lies below the line through the end points $(a, f (a))$ and $(b, f (b))$.

For historical reasons, optimization problems are often referred to as programming problems, which are not related to
computer programming.

Optimization of continuous and smooth functions are closely related to nonlinear equation solving,
because extremal values of a function $f (x)$ correspond to points where its derivative, or gradient, is zero.

Finding candidates for the optimal value of $f(x)$ is therefore equivalent to solving the (in general nonlinear)
equation system $\nabla f (x) = 0$. 

However, a solution to $\nabla f (x) = 0$, which is known as a stationary point, does not necessarily correspond to a minimum of $f (x)$; it can also be maximum or a saddle point.

Candidates obtained by solving $\nabla f (x) = 0$ should therefore be tested for optimality. For unconstrained objective functions the higher-order derivatives, or Hessian matrix

$$ {H_f (x)}_{ij} = \frac{\partial^2 f(x)}{\partial x_i \partial x_j} $$

for the multivariate case, can be used to determine if a stationary point is a local minimum or not.


- In particular if the second-order derivative is positive, or the Hessian positive definite, when evaluated at stationary point $x^*$, then $x^*$ is a local minimum. 

- Negative second-order derivative, or negative definite
Hessian, correspond to a local maximum and 

- A zero second-order derivative, or an indefinite Hessian,
correspond to saddle point.



Algebraically solving the equation system $\nabla f (x) = 0$ and test the candidate solutions for optimality is
therefore one possible strategy for solving an optimization problem. However, it is not always a feasible
method. In particular, we may not have an analytical expression for $f (x)$ from which we can compute the
derivatives, and the resulting nonlinear equation system may not be easy to solve, especially not to find all of
its roots. For such cases, there are alternative numerical optimization approaches,

## univariate optimization

Problem:

$$f(x) = x^2 + 5  $$
$$\min_x f(x) = ?$$


```python
var("X,Y")
```


```python
Y = X**2 + 5
plot(Y)
```


```python
dY = Y.diff()
pc = solve(dY)
pc
```


```python
ddY = dY.diff()
ddY.subs(X,pc[0])
```

## Newton’s method

$$x_{k+1} = x_k - f'(x_k) / f''(x_k) $$

This formula also requires
evaluating both the derivative and the second-order derivative in each iteration. If analytical expressions for
these derivatives are available, this can be a good method. If only function evaluations are available, the
derivatives may be approximated using an analog of the secant method for root finding.

Problem
$$f(x) = x^2 + 5  $$
$$\min_x f(x) = ?$$


```python

```


```python
f = lambda x: x **2  + 5
df = lambda x: 2*x  
ddf = lambda x: 2

x = np.arange(-5,5,.1)

plt.plot(x,f(x))

xk = -4 
for k in range(10):
    xk = xk - df(xk)/ddf(xk)
xk
```


```python
f = lambda x: x **2 - 5*x + 5
df = lambda x: 2*x -5 
ddf = lambda x: 2

x = np.arange(-5,5,.1)

plt.plot(x,f(x))

xk = -4 
for k in range(10):
    xk = xk - df(xk)/ddf(xk)
xk
```




```python
var("X,Y")
f = X**2 + Y**2
f
```


```python
#gradient
nablaf = [ f.diff(x) for x in [X,Y] ]
nablaf
```


```python
solve(nablaf)
```


```python
#hessian
[[f.diff(x).diff(y) for x in [X,Y] for y in [X,Y]]]
```


```python
H= hessian(f,[X,Y])
H
```




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




```python
H.det()
```


```python
a =H.subs( {X:0,Y:0}).eigenvals()
a
#if H(x^*) is positive definite then x^* is a local minimum
#if H(x^*) is negative definite then x^* is a local maximum
```



## Ejemplo


```python
f = X**3 + Y**3 - 3*X*Y
f
```


```python
#gradient
nablaf = [ f.diff(x) for x in [X,Y] ]
nablaf
```


```python
solve(nablaf)
```


```python
H= hessian(f,[X,Y])
H
```




$$\left[\begin{matrix}6 X & -3\\-3 & 6 Y\end{matrix}\right]$$




```python
a =H.subs({X:0,Y:0}).eigenvals()
a
```


```python
a =H.subs({X:1,Y:1}).eigenvals()
a
```


```python

```


```python

```

## Ejemplo

$f(x_1,x_2) = e^{-(x_1^2 + x_2^2) }$


```python
f = exp(-X**2 - Y**2)
f
```


```python
#gradient
nablaf = [ f.diff(x) for x in [X,Y] ]
solve(nablaf)
```


```python
H= hessian(f,[X,Y])
H
```




$$\left[\begin{matrix}4 X^{2} e^{- X^{2} - Y^{2}} - 2 e^{- X^{2} - Y^{2}} & 4 X Y e^{- X^{2} - Y^{2}}\\4 X Y e^{- X^{2} - Y^{2}} & 4 Y^{2} e^{- X^{2} - Y^{2}} - 2 e^{- X^{2} - Y^{2}}\end{matrix}\right]$$




```python
a =H.subs({X:0,Y:0}).eigenvals()
a
```


```python

x = np.arange(-2, 2, 0.25)
y = np.arange(-2, 2, 0.25)

xx, yy = np.meshgrid(x, y)
F = np.exp(-xx**2 - yy**2)

fig = plt.figure()
ax = fig.gca(projection="3d")
surf = ax.plot_surface(xx, yy, F, rstride=1, cstride=1, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False)

ax.set_zlim(-1.01, 1.01)
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.show()
```

$$f(x_1,x_2) = x_1^2 - x_2^2 $$


```python
f = X**2 - Y**2
f
```


```python
#gradient
nablaf = [ f.diff(x) for x in [X,Y] ]
solve(nablaf)
```


```python
H= hessian(f,[X,Y])
H
```




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




```python
a =H.subs({X:0,Y:0}).eigenvals()
a
```


```python

```

## Ejemplo

$$f(x_1,x_2) = -x_1x_2 e^{-(x_1^2 + x_2^2)/2 }$$


```python
f = -X*Y*exp( -(X**2 + Y **2)/2 )
f
```


```python
#gradient
nablaf = [ f.diff(x) for x in [X,Y] ]
solve(nablaf)
```


```python
H= hessian(f,[X,Y])
H
```




$$\left[\begin{matrix}- X^{3} Y e^{- \frac{X^{2}}{2} - \frac{Y^{2}}{2}} + 3 X Y e^{- \frac{X^{2}}{2} - \frac{Y^{2}}{2}} & - X^{2} Y^{2} e^{- \frac{X^{2}}{2} - \frac{Y^{2}}{2}} + X^{2} e^{- \frac{X^{2}}{2} - \frac{Y^{2}}{2}} + Y^{2} e^{- \frac{X^{2}}{2} - \frac{Y^{2}}{2}} - e^{- \frac{X^{2}}{2} - \frac{Y^{2}}{2}}\\- X^{2} Y^{2} e^{- \frac{X^{2}}{2} - \frac{Y^{2}}{2}} + X^{2} e^{- \frac{X^{2}}{2} - \frac{Y^{2}}{2}} + Y^{2} e^{- \frac{X^{2}}{2} - \frac{Y^{2}}{2}} - e^{- \frac{X^{2}}{2} - \frac{Y^{2}}{2}} & - X Y^{3} e^{- \frac{X^{2}}{2} - \frac{Y^{2}}{2}} + 3 X Y e^{- \frac{X^{2}}{2} - \frac{Y^{2}}{2}}\end{matrix}\right]$$




```python
#H.subs({X:0,Y:0}).eigenvals() #silla
#H.subs({X:-1,Y:-1}).eigenvals() #Min
H.subs({X:1,Y:1}).eigenvals() #Min
#H.subs({X:-1,Y:1}).eigenvals() #Max
#H.subs({X:1,Y:-1}).eigenvals() #Max
```


```python

```


```python
%matplotlib notebook
```


```python
x = np.arange(-3, 3, .1)
y = np.arange(-3, 3, .1)

xx, yy = np.meshgrid(x, y)
F = -xx*yy*np.exp( -(xx**2 + yy**2)/2 )

fig = plt.figure()
ax = fig.gca(projection="3d")
surf = ax.plot_surface(xx, yy, F, rstride=1, cstride=1, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False)

ax.set_zlim(-1.01, 1.01)
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.show()
```


    <IPython.core.display.Javascript object>







```python

```


```python
f = (X -1)**4 + 5*(Y -1)**2 - 2*X*Y
#gradient
nablaf = [ f.diff(x) for x in [X,Y] ]
sol =solve(nablaf)
sol
```


```python
_[0]
```


```python
H= hessian(f,[X,Y])
H
```




$$\left[\begin{matrix}12 \left(X - 1\right)^{2} & -2\\-2 & 10\end{matrix}\right]$$




```python
H.subs({X:1,Y:1}).eigenvals() #Min
```


```python

```


```python

```

http://www.dtic.upf.edu/~gharo/anum/apunts/hessian.pdf
http://mat.gsia.cmu.edu/classes/QUANT/NOTES/chap3.pdf

http://www2.econ.iastate.edu/classes/econ500/hallam/documents/Opt_Simple_Multi_000.pdf

e3a349ccf7f77a8edb5e30c6fc653fea470edeb4
ipynb
Jupyter Notebook
Modulo3Simat/3_MaxMinSympy/.ipynb_checkpoints/Optimizacion-checkpoint.ipynb
IntroCursos/MachingLearning
a3d858c8c8fbb8b6d3a1cb445cf5045e9971930c
MIT
Modulo3Simat/3_MaxMinSympy/.ipynb_checkpoints/Optimizacion-checkpoint.ipynb
IntroCursos/MachingLearning
a3d858c8c8fbb8b6d3a1cb445cf5045e9971930c
MIT
Modulo3Simat/3_MaxMinSympy/.ipynb_checkpoints/Optimizacion-checkpoint.ipynb
IntroCursos/MachingLearning
a3d858c8c8fbb8b6d3a1cb445cf5045e9971930c
MIT
Qwen/Qwen-72B
1. YES
2. YES