<a href="https://colab.research.google.com/github/dalexa10/EngineeringDesignOptimization/blob/main/EngineeringOptimization/Local_Implementation/lab02.ipynb" target="_parent"></a>

<div>

</div>
<center><h1>Mathematical Optimization for Engineers</h1></center>
<center><h2>Lab 2</h2></center>
<center><h2>Basic math</h2></center>

$$\newcommand{\mr}[1]{\mathrm{#1}}
\newcommand{\D}{\displaystyle}
\newcommand{\bm}[1]{\text{\mathbf $#1$}}
\newcommand{\bx}{\mathbf{ x}}
\newcommand{\f}{\mathbf{{f}}}
\newcommand{\g}{\mathbf{ g}}
\newcommand{\h}{\mathbf{ h}}
\newcommand{\R}{\mathbb R}
\newcommand{\A}{\mathbf{ A}}
\newcommand{\br}{\boldsymbol{r}}
\newcommand{\bp}{\boldsymbol{p}}
\newcommand{\bnabla}{\mathbf{\nabla}}
$$
In this lab, we will learn about Jacobians, gradients and Hessian matrices.

<u>Notation</u>: Please note that throughout the course, we will denote matrices and vectors with boldface letters. Their components will be denoted by normal letters with subscripts. For example,
$$
\bx = \left(\begin{array}{c}
x_1 \\ 
\vdots \\ 
x_n
\end{array} \right)
$$

## Jacobian

Let $\ \mathbf \f:\R^n \rightarrow \R^m,\,\bx \mapsto \f(\bx)$ be a continuously differentiable function, where
$$
\bx = \left(\begin{array}{c}
x_1 \\ 
\vdots \\ 
x_n
\end{array} \right) , \qquad \f(\bx)  = \left(\begin{array}{c}
f_1(\bx) \\ 
\vdots \\ 
f_m(\bx) 
\end{array} \right).
$$
The Jacobian $\f'(\bx)$ is defined by the matrix
$$
\f'(\bx) = \left(  \begin{array}{ccc}
\frac{\partial f_1(\bx)}{\partial x_1} & \cdots & \frac{\partial f_1(\bx)}{\partial x_n} \\ 
\vdots & \ddots & \vdots \\ 
\frac{\partial f_m(\bx)}{\partial x_1} & \cdots & \frac{\partial f_m(\bx)}{\partial x_n}
\end{array}  \right).
$$

## Gradient

Now, let $f:\R^n \rightarrow \R,\,\bx \mapsto f(\bx)$ be a scalar-valued continuously differentiable function
with vector-valued arguments.
The gradient $\bnabla f(\bx)$ and Jacobian $f'(\bx)$ are defined as the column vector and row vector, respectively,

$$
\bnabla f(\bx)  = \left(\begin{array}{c}
	\frac{\partial f(\bx)}{\partial x_1} \\ 
	\vdots \\ 
	\frac{\partial f(\bx)}{\partial x_n}
\end{array} \right), \qquad  f'(\bx)  = \left(\begin{array}{ccc}  
\frac{\partial f(\bx)}{\partial x_1} & 
\cdots & 
\frac{\partial f(\bx)}{\partial x_n}
\end{array}\right).
$$

Note that  $\bnabla f(\bx)^T=f'(\bx)$.<br>
<br>
This is also valid for any general $\ \mathbf \f:\R^n \rightarrow \R^m,\,\bx \mapsto \f(\bx)$

### Gradient of the scalar product of two functions

Let $\g,\h:\R^n\rightarrow\R^n$ be two continuously differentiable functions. We want to compute the
gradient of $f:\R^n\rightarrow \R,\;\bx \mapsto f(\bx) = \g(\bx)^T\h(\bx)$, where $\bx\in\R^n$.
We have

$$
f(\bx) = \g(\bx)^T\h(\bx) = \sum_{i=1}^{n} g_i(\bx)h_i(\bx).
$$

The derivative with respect to $x_j$ ($1 \le j \le n$) is computed by the application of the product rule

\begin{equation}\label{eq:1}
\frac{\partial f(\bx)}{\partial x_j} = \sum_{i=1}^{n} \left(\frac{\partial g_i(\bx)}{\partial x_j}h_i(\bx)  +g_i(\bx)\frac{\partial h_i(\bx)}{\partial x_j}\right).
\end{equation}

With the notations

$$
\frac{\partial \g(\bx)}{\partial x_j}   = \left(\begin{array}{c}
\frac{\partial g_1(\bx)}{\partial x_j} \\ 
\vdots \\ 
\frac{\partial g_n(\bx)}{\partial x_j}
\end{array} \right) \quad \text{and} \quad 
\frac{\partial \h(\bx)}{\partial x_j}   
= \left(\begin{array}{c}
\frac{\partial h_1(\bx)}{\partial x_j} \\ 
\vdots \\ 
\frac{\partial h_n(\bx)}{\partial x_j}
\end{array} \right), \text{ respectively},
$$

we can rewrite the equation as

$$
\frac{\partial f(\bx)}{\partial x_j} = \frac{\partial \g(\bx)}{\partial x_j} ^T \h(\bx) + 
\g(\bx)^T \frac{\partial \h(\bx)}{\partial x_j}  = \h(\bx)^T\frac{\partial \g(\bx)}{\partial x_j}
+\g(\bx)^T \frac{\partial \h(\bx)}{\partial x_j}
$$

Finally,

$$
\bnabla f(\bx)^T = f'(\bx) = \h(\bx)^T \g'(\bx) + \g(\bx)^T\h'(\bx).
$$

$$
\implies \bnabla f(\bx) = \bnabla g(\bx) \ \h(\bx) + \bnabla h(\bx) \ \g(\bx).
$$

### Derivative of quadratic form
If $\A\in\R^{n\times n}$ is a symmetric matrix, the function

$$
f:\R^n\rightarrow \R,\quad\bx\mapsto f(\bx) = \bx^T\,\A\,\bx
$$ is called a quadratic form. <br>
<br>
With the definitions,

$$
\g(\bx) := \bx \ \text{and } \ \h(\bx):= \A\,\bx,
$$

we have $f(\bx) = \g(\bx)^T\h(\bx)$, i.e., exactly the situation as above.
With $\g'(\bx) = \mathbf{ I}$, where $\mathbf{ I}$ denotes the unity matrix, and $\h'(\bx) = \A$, it is more or less easy to see that the gradient of $f$ is given by

$$
\bnabla f(\bx) = 2\,\A\,\bx.
$$

### Example 1
Let the function $f:\R^2 \rightarrow \R$ be defined as
$ f(x,y) = \frac{1}{2} (x^2 + \alpha y^2), \alpha \in \R $

Find all stationary points of the function $f$.
<details> 
We compute the gradient and set it equal to zero:
$
\begin{equation*}
\frac{\partial f}{\partial x} = x, \; \; \frac{\partial f}{\partial y} =
\alpha y.
\end{equation*}
$
   
If $\alpha \neq 0$ then the stationary point is $x^{*}=(0,0)$.	
If $\alpha = 0$ then any $y$ satisfies $\frac{\partial f}{\partial y} =
0$. Hence we have a "line" of stationary points: $(0,u), \; u \in
\mathbb{R}$.
</details>

Calculate the Hessian of the function $f$ for arbitrary $x$ and $y$.
<details>
We have a smooth function (at least twice continuous and differentiable). Thus, we compute:
$\begin{equation*}
	\frac{\partial^2 f}{\partial x^2} = 1, \; \; \frac{\partial^2 f}{\partial
		y^2} = \alpha, \;\;  \frac{\partial^2 f}{\partial x \partial y} = 0.
\end{equation*} $
    
to obtain the Hessian:
$\begin{equation}
	H := \begin{pmatrix}
	1  & 0 \\
	0 & \alpha
	\end{pmatrix}.
	\end{equation}$
</details>

What are the eigenvalues of the Hessian of the function $f$ with respect to $x$ and $y$?
<details>
  Since the Hessian is a diagonal matrix one can easily retrieve
	its eigenvalues as the elements on the diagonal: $\{1, \alpha\}$.
</details>

Characterize the stationary points for positive and negative $\alpha$.
<details>
If $\alpha > 0$ then $H$ is positive definite, hence by means of the second order sufficient conditions we conclude that
	$(0,0)$ is a minimum point.
	
If $\alpha < 0$ then $H$ is indefinite, rendering us the saddle
	point $(0,0)$.
</details>

Characterize the convexity of the function for every $\alpha$.
<details>
The convexity can be determined by means of the definiteness of the Hessian $H$:
    
a. for $\alpha > 0$ $H$ is positive definite $\Rightarrow$ $f$ is strictly convex;
    
b. for $\alpha =0$ $H$ is positive semidefinite	$\Rightarrow$ $f$ is convex (but not strictly);
    
c. for $\alpha <0$ $H$ is indefinite 
		$\Rightarrow$ $f$ is neither convex nor concave.    
</details>


### Example 2: Rosenbrock function
The Rosenbrock function is a famous test function used in optimization. <br>
<br> 
It is defined by: 
$f:\R^2 \rightarrow \R$
$ f(x,y) = (a-x)^2 + b(y-x^2)^2 $ with $a=1, b=100$

Find all stationary points of the function $f$.
<details> 
We need to set the gradient to zero. 
    
$\begin{align}
\frac{\partial f}{\partial x} =&\,\, 2(-1+x) + 200(y-x^2)(-2x) = 0,\label{eq:rosenbgradx} \\
\frac{\partial f}{\partial y} =&\,\,  200(y-x^2) = 0. \label{eq:rosenbgrady}
\end{align}
$
    
Substituting the second equation into the first, we obtain: 
$\begin{align}
2x-2 = 0
\end{align}$
It follows: 
$\begin{align}
x^*=1\\
y^*=1
\end{align}$
    
    
</details>

Calculate the Hessian of the function $f$ for the stationary points found. (Hint: use Python)
<details> 
We compute the Hessian as: 
$\begin{pmatrix}
400(3x^2-y)+2 & -400x\\
-400x & 200
\end{pmatrix}
$


For our stationary point, we have 
$\begin{pmatrix}
802 & -400\\
-400 & 200
\end{pmatrix}
$.



</details>

What are the eigenvalues of the Hessian of the function $f$? (Hint: use Python)
<details> 
The eigenvalues are 0.4 and 1000. Hence, the Hessian is positive definite. 
The second order sufficient conditions for a minimum are fulfilled. 
   
</details>



```python
def rosenbrock(x):
    return ((x[0]-1)**2 + 100*(x[1]-x[0]**2)**2)
```


```python
# compute hessian using autograd
import autograd
import numpy as np
x0 = np.array([1.,1.])
hessian_rosenbrock = autograd.hessian(rosenbrock)(x0)
print(hessian_rosenbrock)
```

    [[ 802. -400.]
     [-400.  200.]]



```python
# compute the eigenvalues using numpy
import numpy as np
w, v = np.linalg.eig(hessian_rosenbrock)
print(w)
```

    [1.00160064e+03 3.99360767e-01]



```python

```

e3c3b6dad5761c68ba099d43a871c73b62991c1f
ipynb
Jupyter Notebook
EngineeringOptimization/Local_Implementation/lab02.ipynb
dalexa10/EngineeringDesignOptimization
eb5b5e4edd773aef629f59aea8a9771af41bd224
MIT
EngineeringOptimization/Local_Implementation/lab02.ipynb
dalexa10/EngineeringDesignOptimization
eb5b5e4edd773aef629f59aea8a9771af41bd224
MIT
EngineeringOptimization/Local_Implementation/lab02.ipynb
dalexa10/EngineeringDesignOptimization
eb5b5e4edd773aef629f59aea8a9771af41bd224
MIT
Qwen/Qwen-72B
1. YES
2. YES