nonlinear_programming.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "366edae6-c36d-4b71-9d27-83de8f4d7ee7",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Nonlinear Programming with Constraints\n",
    "\n",
    "## Anselm Hudde\n",
    "\n",
    "![](optimization.png)\n",
    "**Source: Venter, G. (2010). Review of optimization techniques.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bda742a-3bc3-4d75-9c5d-867099ff875c",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%run nonlinear_programming.py"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25e78631",
   "metadata": {},
   "source": [
    "## What is nonlinear programming about?\n",
    "\n",
    "The goal of optimization is to find the minimum (or maximum, resp.) value of a function $f$.\n",
    "\n",
    "- How does a company maximize its revenue, considering the available resources?\n",
    "\n",
    "- Given specific return goals, what asset allocation minimizes the risk?\n",
    "\n",
    "- How to calibrate model parameters (statistical models, processes to model shares or interest rates, neural networks)?\n",
    "\n",
    "## Asset allocation example\n",
    "\n",
    "An investor wants to invest € $10\\,000$, with a portfolio which consists of two shares, $A$ and $B$, respectively.\n",
    "\n",
    "- Expected yearly returns:\n",
    "\n",
    "    - Share $A$: $r_{A} = 7\\%$\n",
    "\n",
    "    - Share $B$: $r_{B} = 9\\%$\n",
    "    \n",
    "    \n",
    "- Volatilities:\n",
    "\n",
    "    - Share $A$: $\\sigma_{A} = 20\\%$\n",
    "\n",
    "    - Share $B$: $\\sigma_{B} = 30\\%$\n",
    "    \n",
    "    \n",
    "- Correlation:\n",
    "\n",
    "    - $\\rho_{A, B} = 0.7$\n",
    "\n",
    "The investor wants to invest her money in such a way, that the expected return is $5\\%$, and where risk (volatility) is minimized.\n",
    "A portfolio consisting of $x_{1} \\times 10\\,000$ € in Share $A$ und $x_{2} \\times 10\\,000$ € in share $B$ can be expressed as\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "f (x_{1}, x_{2})\n",
    "=\n",
    "\\sqrt{ \\sigma_{A}^2 x_1^2 + \\sigma_{B}^2 x_2^2 + 2 \\sigma_{A} x_1 \\sigma_{B} x_2 \\rho_{A, B}}.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Our goal is\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\text{Minimize }& f(x) & ~~~~~~~~ &\\text{ (the volatility)}\\\\\n",
    "&\\text{where }& x_{1} + x_{2} - 1 \\leq 0 & &\\text{ (do not invest more than 100 %)}\\\\\n",
    "&\\text{and }& r_{A} x_{1} + r_{B} x_{2} - 0.05 = 0 & &\\text{ (5 % expected return)}.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "By setting\n",
    "\n",
    "\\begin{equation}\n",
    "g(x_1, x_2) = x_1 + x_2 - 1\n",
    "\\end{equation}\n",
    "and\n",
    "\\begin{equation}\n",
    "h(x_1, x_2) = r_{A} x_{1} + r_{B} x_{2} - 0.05,\n",
    "\\end{equation}\n",
    "\n",
    "we can regard our problem as instance of a more general formalization of a nonlinear optimization problem with constraints.\n",
    "\n",
    "## Nonlinear programming with constraints\n",
    "\n",
    "Let $D \\subseteq \\mathbb{R}^n$ be a subset, and let\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "f \\colon D \\to \\mathbb{R} \\\\\n",
    "g \\colon D \\to \\mathbb{R} \\\\\n",
    "h \\colon D \\to \\mathbb{R}\n",
    "\\end{split}\n",
    "\\end{equation}\n",
    "\n",
    "be three (possibly nonlinear) functions.\n",
    "\n",
    "**Goal of nonlinear programming:**\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\text{Minimize }& f(x) \\\\\n",
    "\\text{where }& g(x) \\leq 0 \\\\\n",
    "\\text{and }& h(x) = 0 \\text{ holds true.}\n",
    "\\end{split}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29bb3a6e-f567-4a67-8a00-e22646a019bb",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "source": [
    "## Example function\n",
    "\n",
    "We consider the function\n",
    "\n",
    "\\begin{equation}f(x_{1}, x_{2}) = \\sin(x_{1}) + 0.05 x_{1}^2 + 0.1 x_{2}^2, \\end{equation}\n",
    "\n",
    "which we want to minimize:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be85e606",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def f(x):\n",
    "    return np.sin(x[0]) + 0.05*x[0]**2 + 0.1*x[1]**2\n",
    "\n",
    "surface_plot = plot()\n",
    "surface_plot.plot_surface(-5,1,-3,3,f)\n",
    "\n",
    "contour_plot = plot()\n",
    "contour_plot.plot_contour(-5,1,-3,3,f)\n",
    "\n",
    "show_plot(contour_plot, surface_plot)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e406e17b-70d2-4f3f-81b2-3354f42acf0e",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## Gradient descent: the idea\n",
    "\n",
    "The general idea behind gradient descent is, to always go in the direction with the steeped ascent (resp. descent).\n",
    "\n",
    "![](optimization.png)\n",
    "\n",
    "But how do we find this direction? To this end, we need to be capable of measuring the slope.\n",
    "\n",
    "This can the done with the **Gradient**:\n",
    "\n",
    "**Definition:**\n",
    "The gradient of the function $f \\colon \\mathbb{R}^{n} \\to \\mathbb{R}$ is the function\n",
    "\n",
    "\\begin{equation}\n",
    "\\nabla f \\colon \\mathbb{R}^{n} \\to \\mathbb{R}^{n};~~\n",
    "x =\n",
    "\\begin{pmatrix}\n",
    "x_{1} \\\\ \\vdots \\\\ x_{n}\n",
    "\\end{pmatrix}\n",
    "\\mapsto\n",
    "\\begin{pmatrix}\n",
    "\\tfrac{\\partial f}{\\partial x_{1}} (x) \\\\ \\vdots \\\\ \\tfrac{\\partial f}{\\partial x_{n}} (x)\n",
    "\\end{pmatrix}.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73cf3abe-4911-45c0-99b7-78eb125d5cce",
   "metadata": {
    "tags": []
   },
   "source": [
    "## **Intuition:** The Gradient always points in the direction of the steepest ascent.\n",
    "\n",
    "**Example 1:** $f(x_{1}, x_{2}) = x_1$, $\\nabla f(x_{1}, x_{2}) = \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7392913a-0e0c-4456-8938-893f47d3cd1d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "surface_plot = plot()\n",
    "surface_plot.plot_surface(0,4,0,4,function = lambda x: x[0])\n",
    "surface_plot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13632e34-1665-4598-8873-2bc804728365",
   "metadata": {},
   "source": [
    "**Example 2:** $f(x_{1}, x_{2}) = \\tfrac{1}{2} (x_{1} + x_{2})$, $\\nabla f(x_{1}, x_{2}) = \\begin{pmatrix} \\tfrac{1}{2} \\\\ \\tfrac{1}{2} \\end{pmatrix}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0cac26d-b0e0-45c8-aa7e-34766f499044",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "surface_plot = plot()\n",
    "surface_plot.plot_surface(0,4,0,4,function = lambda x: (1/3)*(x[0]+x[1]))\n",
    "surface_plot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3635e17b-370a-41e8-9254-f9d3139664f6",
   "metadata": {},
   "source": [
    "## Calculating the gradient of $f$\n",
    "\n",
    "The gradient of\n",
    "$𝑓(x_{1}, x_{2})= \\sin(𝑥_{1})+0.05 𝑥_{1}^{2} + 0.1 x_{2}^{2}$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\nabla f(x_{1}, x_{2})\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "\\frac{\\partial f}{\\partial x_{1}} \\\\ \\frac{\\partial f}{\\partial x_{2}} \\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "\\cos(x_{1}) + 0.1 x_{1} \\\\ 0.2 x_{2}\n",
    "\\end{pmatrix}\n",
    "\\end{equation}\n",
    "\n",
    "**Examples:**\n",
    "\n",
    "\\begin{equation}\n",
    "\\nabla f(-4, -2) =\n",
    "\\begin{pmatrix}\n",
    "-1.05 \\\\ -0.40\n",
    "\\end{pmatrix},\n",
    "\\end{equation}\n",
    "\\begin{equation}\n",
    "\\nabla f(-1, 2) =\n",
    "\\begin{pmatrix}\n",
    "0.44 \\\\ 0.20\n",
    "\\end{pmatrix}.\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82db4c94-24fe-4804-90df-351454a787df",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def grad_f(x):\n",
    "    return np.array([np.cos(x[0])+0.1*x[0], 0.2*x[1]])\n",
    "\n",
    "contour_plot.add_gradients(grad_f)\n",
    "contour_plot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5117f6f",
   "metadata": {},
   "source": [
    "## The steps of the gradient descent\n",
    "\n",
    "We start with a point $x^{(0)} \\in \\mathbb{R}^{2}$.\n",
    "For the $n$-th step, we calculate the gradient\n",
    "$\\nabla f(x^{(n)})$,\n",
    "and go one step in the opposite direction of the gradient, multiplied by the learning rate $\\gamma$:\n",
    "\n",
    "\\begin{equation}\n",
    "\tx^{(n+1)}= x^{(n)} - \\gamma \\nabla f({x}^{(n)}),\\ n\\geq 0.\n",
    "\\end{equation}\n",
    "\n",
    "### Learning rate $\\gamma = 1$:\n",
    "\n",
    "We start with the point\n",
    "$x^{(0)} = \\begin{pmatrix} -4 \\\\ -2 \\end{pmatrix}$\n",
    "and apply a learning rate of $\\gamma = 1$.\n",
    "Then we have\n",
    "\n",
    "\\begin{equation} \\nabla f(x^{(0)}) =\n",
    "\\begin{pmatrix}\\cos(-4)+0.1\\cdot(-4) \\\\ 0.2\\cdot(-2) \\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix} -1.05 \\\\ -0.40 \\end{pmatrix}\n",
    "\\end{equation}\n",
    "\n",
    "and therefore\n",
    "\n",
    "\\begin{equation}\n",
    "x^{(1)} = x^{(0)} + \\gamma_{0} \\cdot \\nabla f(x^{(0)})\n",
    "=\n",
    "\\begin{pmatrix} -4 \\\\ -2 \\end{pmatrix}\n",
    "- 1 \\cdot \n",
    "\\begin{pmatrix} -1.05 \\\\ -0.40 \\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "-2.95 \\\\ -1.60\n",
    "\\end{pmatrix}.\n",
    "\\end{equation}\n",
    "The next steps follow analogously:\n",
    "\n",
    "$x^{(2)} = \\begin{pmatrix} -1.67 \\\\ -1.28\\end{pmatrix}, ~\n",
    "x^{(3)} = \\begin{pmatrix} -1.40 \\\\ -1.02\\end{pmatrix}, ~\\dots$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ea11cba-4570-4aaf-b9cf-6e1f4b898999",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Gradient descent with learning rate gamma = 1:\n",
    "\n",
    "def f(x):\n",
    "    return np.sin(x[0]) + 0.05*x[0]**2 + 0.1*x[1]**2\n",
    "\n",
    "def grad_f(x):\n",
    "    return np.array([np.cos(x[0])+0.1*x[0], 0.2*x[1]])\n",
    "\n",
    "contour_plot = plot()\n",
    "contour_plot.plot_contour(-5,1,-3,3,f)\n",
    "\n",
    "surface_plot = plot()\n",
    "surface_plot.plot_surface(-5,1,-3,3,f, opacity=0.5)\n",
    "\n",
    "i = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "041ea80e",
   "metadata": {},
   "source": [
    "The following cell executes the single steps of gradient descent with the learning rate $\\gamma = 1$.\n",
    "Running this cell repeatedly shows, how the gradient descent converges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a5617af",
   "metadata": {},
   "outputs": [],
   "source": [
    "contour_plot.add_gradient_descent(x0=[-4, -2], function = f, grad=grad_f, gamma=1, iterations=i, color = \"#636EFA\")\n",
    "surface_plot.add_gradient_descent_surface(x0=[-4, -2], function = f, grad=grad_f, gamma=1, iterations=i, color = \"#636EFA\")\n",
    "show_plot(contour_plot, surface_plot)\n",
    "i += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d308b5e-5660-44a2-aabe-cf2e385b9d8a",
   "metadata": {},
   "source": [
    "### Gradient descent with a learning rate of $\\gamma = 0.1$ is slower:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "031e363b-6d69-4bd9-84a6-9ca6a5d5fed8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "i = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d16681b",
   "metadata": {},
   "source": [
    "The following cell executes the single steps of gradient descent with the learning rate $\\gamma = 0.1$.\n",
    "Running this cell repeatedly shows, how the gradient descent converges slowly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "669954b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gradient descent with learning rate gamma = 0.1:\n",
    "\n",
    "contour_plot.add_gradient_descent(x0=[-4, -2], function = f, grad=grad_f, gamma=0.1, iterations=i, color = \"#EF553B\")\n",
    "surface_plot.add_gradient_descent_surface(x0=[-4, -2], function = f, grad=grad_f, gamma=0.1, iterations=i, color = \"#EF553B\")\n",
    "show_plot(contour_plot, surface_plot)\n",
    "i+=1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aba80109",
   "metadata": {},
   "source": [
    "### Gradient descent with a learning rate of $\\gamma = 2$ does not converge:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7bd4ddd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "i = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa58aa31",
   "metadata": {},
   "source": [
    "The following cell executes the single steps of gradient descent with the learning rate $\\gamma = 0.1$.\n",
    "Running this cell repeatedly shows, that is this case, gradient descent does not converge.\n",
    "The solution jumps back and forth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62039c75",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gradient descent with learning rate gamma = 2:\n",
    "\n",
    "contour_plot.add_gradient_descent(x0=[-4, -2], function = f, grad=grad_f, gamma=2, iterations=i, color = \"#00CC96\")\n",
    "surface_plot.add_gradient_descent_surface(x0=[-4, -2], function = f, grad=grad_f, gamma=2, iterations=i, color = \"#00CC96\")\n",
    "show_plot(contour_plot, surface_plot)\n",
    "i += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a48f4425-5044-4b9e-81f2-70385c1a0eca",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "source": [
    "## Impact of the initial guess:\n",
    "\n",
    "If the target function has several local minima, the choice of the initial guess also plays a vital role:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1ff7250-8aed-4c9b-a5ee-1447587bbd65",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "def f(x):\n",
    "    return np.sin(x[0]) + 0.05*x[0]**2 + 0.1*x[1]**2\n",
    "\n",
    "contour_plot = plot()\n",
    "contour_plot.plot_contour(-5,6,-3,3,f)\n",
    "\n",
    "surface_plot = plot()\n",
    "surface_plot.plot_surface(-5,6,-3,3,f)\n",
    "\n",
    "show_plot(contour_plot, surface_plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e5647d5",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "x0=[random.uniform(-5,6),random.uniform(-3,3)]\n",
    "\n",
    "contour_plot.add_gradient_descent(x0=x0, function = f, grad=grad_f, gamma=1, iterations=30)\n",
    "surface_plot.add_gradient_descent_surface(x0=x0, function = f, grad=grad_f, gamma=1, iterations=30)\n",
    "show_plot(contour_plot, surface_plot)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39148f36-835e-4ef2-9073-5764c20f03f1",
   "metadata": {},
   "source": [
    "## Nonlinear programming with `SciPy`\n",
    "\n",
    "The presented method is still too simple for practical use, in practical application there are still many problems to be considered.\n",
    "Therefore, one usually takes an already existing implementation.\n",
    "The Python package `SciPy` with the function `minimize` is very suitable for this.\n",
    "It suffices to enter the function to be minimized and the initial guess:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1739a455",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Optimizing with SciPy\n",
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "def f(x):\n",
    "    return np.sin(x[0]) + 0.05*x[0]**2 + 0.1*x[1]**2\n",
    "\n",
    "minimize(fun=f, x0=np.array([-4,-2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d01f650-5cc6-4e3f-aca4-f3e8f6f5babd",
   "metadata": {},
   "source": [
    "## Constraints of the form $g(x) \\leq 0$\n",
    "\n",
    "An example for a constraint of the form $g(x) \\leq 0$ is\n",
    "\n",
    "\\begin{equation}\n",
    "-(x_{1}+2)^{2} + x_{2}^{3} \\leq 0.\n",
    "\\end{equation}\n",
    "\n",
    "This inequality defines a domain, in which the solution lives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2160642d-9121-4f7c-86dc-9ebc109f8efd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def h(x):\n",
    "    return (x[0]+3)**3-x[1]\n",
    "\n",
    "def g_(x):\n",
    "    return - ((x+2)**2)**(1/3)\n",
    "\n",
    "x1 = np.linspace(-5, -2, 100).tolist()\n",
    "y1 = [np.sqrt(np.abs(2.25 - (x+3.5)**2)) for x in x1]\n",
    "x2 = np.linspace(-2, -5, 100).tolist()\n",
    "y2 = [-np.sqrt(np.abs(2.25 - (x+3.5)**2)) for x in x2]\n",
    "x = x1 + x2\n",
    "y = y1 + y2\n",
    "\n",
    "contour_plot = plot()\n",
    "contour_plot.plot_contour(-5,1,-3,3,f)\n",
    "\n",
    "contour_plot.add_trace(go.Scatter(x=x, y=y, fill='tonexty', showlegend = False, marker = {'color' : '#19D3F3'})).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29e55d28-a0cf-427a-a31e-b2c98f46150c",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Constraints of the form $h(x) = 0$\n",
    "\n",
    "A constraint of the form $h(x) = 0$ usually defines a line on which the solution lives.\n",
    "This is usually a restrictive condition, i.e., a worse result is obtained than without a constraint.\n",
    "We show as an example a constraint of the form\n",
    "\n",
    "\\begin{equation}\n",
    "h(x_{1}, x_{2}) = (x_{1} + 3)^2 - x_{2}.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f321586",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "contour_plot.add_h()\n",
    "contour_plot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdc6e338-7adf-4cd4-b108-6e996c3d84bc",
   "metadata": {},
   "source": [
    "## The quadratic penalty method\n",
    "\n",
    "We know now, how to solve unconstrained nonlinear programming optimization problems.\n",
    "But how can we consider constraints, i.e.,\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\text{Minimize }& f(x) \\\\\n",
    "\\textbf{where }& \\mathbf{g(x) \\leq 0} \\\\\n",
    "\\textbf{and }& \\mathbf{h(x) = 0} \\textbf{ holds true?}\n",
    "\\end{split}\n",
    "\\end{equation}\n",
    "\n",
    "One possibility is to add a penalty term to the target function $f$.\n",
    "This penalty term penalizes deviations from the constraints.\n",
    "Such methods are called **penalty methods**.\n",
    "Here, we introduce the **quadratic penalty method**.\n",
    "The problem can be formulated as:\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{Minimize } f(x) + \\alpha \\max(0, g(x))^{2} + \\alpha h(x)^{2}, ~ \\alpha > 0.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f241e9da-308c-4e10-9e7b-02727bead3f0",
   "metadata": {
    "tags": []
   },
   "source": [
    "## The penalty method for $h(x) = 0$:\n",
    "\n",
    "Reminder:\n",
    "\n",
    "\\begin{equation}\n",
    "h(x_{1}, x_{2}) = (x_{1} + 3)^2 - x_{2}.\n",
    "\\end{equation}\n",
    "\n",
    "We choose an $\\alpha > 0$.\n",
    "Then we have to minimize the function\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "f_{\\text{pen}} (x_{1}, x_{2})\n",
    "&=\n",
    "f(x_{1}, x_{2}) + \\alpha h(x_{1}, x_{2})^2.\n",
    "\\end{split}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d394b1f5-1c79-42f8-8338-7f6d3f203451",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "alpha = 0.01\n",
    "\n",
    "def h(x):\n",
    "    return (x[0]+3)**2 - x[1]\n",
    "\n",
    "def f_pen(x):\n",
    "    return f(x) + alpha*h(x)**2\n",
    "\n",
    "surface_plot = plot()\n",
    "surface_plot.plot_surface(-5,1,-3,3,f,opacity = 0.8, colorscale = 'blues')\n",
    "surface_plot.plot_surface(-5,1,-3,3,f_pen, opacity = 1, showscale=False, colorscale = 'oranges')\n",
    "surface_plot.show()\n",
    "surface_plot.write_html(\"bla.html\")\n",
    "\n",
    "def grad_h_sq(x):\n",
    "    return np.array([4*x[0]**3 + 36*x[0]**2 + 108*x[0] + 108 - 4*x[0]*x[1] - 12*x[1],\n",
    "                     -2*x[0]**2 - 12*x[0] + 2*x[1] - 18])\n",
    "\n",
    "contour_plot = plot()\n",
    "contour_plot.plot_contour(-5,1,-3,3,f)\n",
    "contour_plot.add_h()\n",
    "\n",
    "contour_plot.result = [0, 3]\n",
    "\n",
    "gamma = 1\n",
    "\n",
    "x_range = 1\n",
    "y_range = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcf10edc",
   "metadata": {},
   "source": [
    "The gradient of $h(x_{1}, x_{2})^2$ is\n",
    "\n",
    "\\begin{equation}\n",
    "\\nabla h^2(x_{1}, x_{2}) =\n",
    "\\begin{pmatrix}\n",
    "4x_{1}^{3} + 36x_{1}^{2} + 108x_{1} + 108 - 4x_{1}x_{2} - 12x_{2} \\\\\n",
    "-2x_{1}^{2} -12x_{1} + 2x_{2} - 18\n",
    "\\end{pmatrix}.\n",
    "\\end{equation}\n",
    "\n",
    "The Gradient of $f_{ \\text{pen}}$ is\n",
    "\n",
    "\\begin{equation}\n",
    "\\nabla f_{ \\text{pen}} = \n",
    "\\nabla f + \\alpha \\nabla h.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17c79ac3",
   "metadata": {},
   "outputs": [],
   "source": [
    "contour_plot.add_gradient_descent(x0=contour_plot.result, function=f, grad=lambda x : (grad_f(x) + alpha*grad_h_sq(x)), gamma=gamma,\n",
    "                                  iterations=100, constraint = h)\n",
    "\n",
    "contour_plot.show()\n",
    "\n",
    "x_result = contour_plot.result[0]\n",
    "y_result = contour_plot.result[1]\n",
    "\n",
    "contour_plot.contour_zoom(x_result-x_range, x_result+x_range, y_result-y_range, y_result+y_range, f)\n",
    "\n",
    "alpha *= 10\n",
    "gamma /= 10\n",
    "\n",
    "x_range /= 1.1\n",
    "y_range /= 1.1\n",
    "\n",
    "x0 = contour_plot.result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e6aca29",
   "metadata": {},
   "source": [
    "## Nonlinear programming with constraints in python\n",
    "\n",
    "Nonlinear programming with constraints also works well in python:\n",
    "\n",
    "### Optimizing with $h(x) = 0$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1194a4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optimizing with SciPy\n",
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "from scipy.optimize import NonlinearConstraint\n",
    "\n",
    "def f(x):\n",
    "    return np.sin(x[0]) + 0.05*x[0]**2 + 0.1*x[1]**2\n",
    "\n",
    "def h(x):\n",
    "    return (x[0]+3)**2 - x[1]\n",
    "\n",
    "constraints = NonlinearConstraint(h, lb = 0, ub = 0)\n",
    "\n",
    "minimize(fun=f, x0=np.array([3,0]), constraints=constraints)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f88e98de",
   "metadata": {},
   "source": [
    "### Optimizing with $h(x) = 0$ and $g(x) \\leq 0$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d472993c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def g(x):\n",
    "    return -(x[0]+2)**2 + x[1]**3\n",
    "\n",
    "constraints = [NonlinearConstraint(h, lb = 0, ub = 0),\n",
    "               NonlinearConstraint(g, lb = -np.inf, ub = 0)]\n",
    "\n",
    "minimize(fun=f, x0=np.array([-4,-2]), constraints=constraints)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bed26dfe",
   "metadata": {},
   "source": [
    "## Outlook: Automatic Differentiation\n",
    "\n",
    "### Numerical computation of gradients:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial f}{\\partial x_{1}} \\approx \\frac{f(x_{1} + h, x_{2}) - f(x_{1}, x_{2})}{h}, ~ \n",
    "\\frac{\\partial f}{\\partial x_{2}} \\approx \\frac{f(x_{1}, x_{2} + h) - f(x_{1}, x_{2})}{h}, ~\n",
    "h\\text{ small}\n",
    "\\end{equation}\n",
    " \n",
    "$\\Rightarrow$ For each partial derivative a separate difference quotient must be formed.\n",
    "\n",
    "- Particularly demanding application of nonlinear optimization: neural networks\n",
    "- Gradient method essential algorithm for training neural networks\n",
    "\n",
    "### Application for neural networks\n",
    "\n",
    "![](simple_fnn_backprop.svg)\n",
    "<figcaption align = \"center\"><b>Graphic shows a neural network together with the loss function. Source: Maren Hackenberg, AG Machine Learning, Institute of Medical Biometry and Statistics (IMBI), University of Freiburg, Germany, Source: </b> <a href=\"https://github.com/maren-ha/Masterarbeit/blob/master/Figures/simple_fnn_backprop.pdf\">github.com/maren-ha/Masterarbeit/blob/master/Figures/simple_fnn_backprop.pdf</a></figcaption>\n",
    "\n",
    "- Loss function $\\mathcal{L}$ depends on many parameters (e.g. $10\\,000 - 10\\,000\\,000$)\n",
    "- numerical computation of the gradient costly and prone to rounding errors\n",
    "- analytical calculation not practicable and dependent on the exact structure of the neural net\n",
    "\n",
    "### Solution: Automatic Differentiation\n",
    "\n",
    "- Automatic differentiation packages can compute the exact gradient of (many) functions written in python\n",
    "- Hence more efficient optimization of high-dimensional problems possible\n",
    "- The combination of gradient descent and automatic differentiation is a central element for the success of many machine learning methods\n",
    "\n",
    "### Example: Increase in efficiency through Automatic Differentiation\n",
    "\n",
    "We optimize the function\n",
    "\n",
    "\\begin{equation}\n",
    "f(x) = \\exp \\left(- \\sum_{i=1}^{n} i \\cdot x_{i} \\right) + \\sum_{i=1}^{n} x_{i}^{2}.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b620aac8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import autograd.numpy as np\n",
    "from autograd import grad\n",
    "from scipy.optimize import minimize\n",
    "from time import time\n",
    "\n",
    "def f(x):\n",
    "    return np.exp(-np.sum(np.array(range(1, len(x)+1) * x))) + np.sum(x**2)\n",
    "\n",
    "x0=np.ones(2)\n",
    "\n",
    "## Without gradient\n",
    "minimize(fun=f, x0=x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4cb0a525",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Calculating the gradient\n",
    "grad_f = grad(f)\n",
    "minimize(fun=f, x0=x0, jac = grad_f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fefea20",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Benchmarks\n",
    "\n",
    "x0=np.ones(200)\n",
    "\n",
    "t0 = time()\n",
    "for i in range(10):\n",
    "    fd = minimize(fun=f, x0=x0)\n",
    "t_fd = (time() - t0) / 10\n",
    "\n",
    "print(\"Function value finite difference: \", fd.fun)\n",
    "print(\"Computing time finite difference: \", t_fd)\n",
    "\n",
    "t0 = time()\n",
    "for i in range(10):\n",
    "    grad_f = grad(f)\n",
    "    autodiff = minimize(fun=f, x0=x0, jac = grad_f)\n",
    "t_autodiff = (time() - t0) / 10\n",
    "\n",
    "print(\"Function value Autodiff: \", autodiff.fun)\n",
    "print(\"Computing time Autodiff: \", t_autodiff)\n",
    "print(\"(Computational time finite difference) / (Computational time Autodiff):\", t_fd / t_autodiff)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abd61764",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "92f23b4ac9a1ac0aabac2a5fe090cf874f268b01de51ad4092840717b54c19ba"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}