{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Import libraries\n", "import numpy as np\n", "from numpy import linalg as LA\n", "\n", "import random\n", "import time\n", "\n", "import matplotlib.pyplot as plt\n", "plt.rcParams['text.usetex'] = True" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## MA934\n", "\n", "## Nonlinear optimisation\n" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "split", "slideshow": { "slide_type": "slide" } }, "source": [ "### Minimisation in 1 dimension\n", "\n", "Given a function, $f(x)$ of a single variable, the task is to find a minimum of $f$. \n", "\n", "An ordered triple of points $(a,c,b)$ is said to *bracket* a minimum of $f$ if $f(c) < f(a)$ and $f(c) < f(b)$.\n", "\n", "Line search: evaluate $f$ at a new point, $x$, to construct a smaller bracketing triple. Iterate until a desired accuracy is reached." ] }, { "cell_type": "markdown", "metadata": { "cell_style": "split" }, "source": [ "\"array\" " ] }, { "cell_type": "markdown", "metadata": { "cell_style": "split" }, "source": [ "Golden section search is a way to organise this search in an optimal way." ] }, { "cell_type": "markdown", "metadata": { "cell_style": "center", "slideshow": { "slide_type": "slide" } }, "source": [ "### Golden section search\n", "\n", "\"array\" \n", "\n", "Most efficient to choose $x$ in the larger of the two subintervals: \n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Golden section search: choosing the new point, $x$" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "split", "slideshow": { "slide_type": "-" } }, "source": [ "**Case 1**: $[a, c]$ is the larger subinterval : new bracketing triple is either $(a,x,c)$ or $(x,c,b)$. Width of the new bracketing triple is independent of which outcome if :\n", "$$\n", "c - a = b - x.\n", "$$\n", "So we choose\n", "$$\n", "x = a + b - c.\n", "$$\n" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "split", "slideshow": { "slide_type": "fragment" } }, "source": [ "**Case 2**: if $[c,b]$ is the larger subinterval: new bracketing triple is either $(a,c,x)$ or $(c,x,b)$. Width of the new bracketing triple is independent of which outcome if :\n", "$$\n", "x - a = b - c.\n", "$$\n", "So we again choose\n", "$$\n", "x = a + b -c.\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "split", "slideshow": { "slide_type": "slide" } }, "source": [ "### Golden section search : choosing initial \"shape\"\n", "\n", "\"array\" \n", "\n", "The idea is to choose $c$ such that the ratio of the width of the shorter subinterval to the width of the longer one remains constant between iterations." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Golden section search : choosing initial \"shape\"\n" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "split" }, "source": [ "Denote:\n", "\n", "* $p$: width of the longer subinterval in the old triple. \n", "* $q$: width of the shorter subinterval in the old triple.\n", "* $r$: width of shorter subinterval in the *new* triple.\n", "\n", "There are several cases:" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "split" }, "source": [ "Case 1:\n", "\n", "- Case 1(A) : $\\frac{r}{p-r} = \\frac{q}{p}$\n", "- Case 1(B) : $\\frac{r}{q} = \\frac{q}{p}$\n", "\n", "Case 2: (get same equations) \n", "- Case 2(A) : $\\frac{r}{q} = \\frac{q}{p}$ \n", "- Case 2(B) : $\\frac{r}{p-r} = \\frac{q}{p}$" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "split" }, "source": [ "\"array\" \n", "\n", "Eliminating $r$ gives\n", "\n", "$$\n", "\\begin{align*}\n", "\\left(\\frac{q}{p} \\right)^2 &+ \\left(\\frac{q}{p} \\right) - 1 = 0 \\\\ \n", "\\Rightarrow \\left(\\frac{q}{p} \\right) &= \\frac{\\sqrt{5}\\pm 1}{2}.\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Golden section search\n", "\n", "We choose the - sign since we assumed $\\frac{q}{p}<1$ ($\\frac{q}{p} = \\phi$, 1/golden ratio). \n", "\n", "Since $x$ is already determined, $x=a+b−c$, if we start with the correct ratio, $\\frac{q}{p}$ , this will be preserved when as we iterate.\n", "\n", "Convergence is exponential in the number of iterations: width of interval after $n$ iterations is $(b-a)\\,\\phi^n$.\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Golden section search: convergence\n", "\n", "Finding the minimum of $f(x)=x^2$ at $x=0$:\n", "\n", "\"array\" " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Golden section search : implementation\n", "\n", "Implementation is simplified using temporary variables $x_1$ and $x_2$ to store the points $x$ and $c$ (in Case 1) or $c$ and $x$ (in Case 2) and $f_1$ and $f_2$ to store the associated function values. \n", "\n", "Regardless of the order in which previous points were evaluated, the new triple will centre on the point with the smallest value of $f$ found so far. Thus by comparing $f_1$ to $f_2$, there are only two cases:\n", "\n", "\n", "\"array\" \n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Minimisation in $\\mathbb{R^n}$ : line minimisation\n", "\n", "Given a function, $f: \\mathbb{R}^n \\to \\mathbb{R}$, a position vector, $\\mathbf{x} \\in \\mathbb{R}^n$, and a direction vector, $\\mathbf{v} \\in \\mathbb{R}^n$, the *line minimiser* of $f$ from $\\mathbf{x}$ in the direction $\\mathbf{v}$ is the point \n", "$$\n", "\\mathbf{x}^* = \\mathbf{x} + \\lambda^*\\,\\mathbf{v}\n", "$$ \n", "where\n", "$$\n", "\\lambda^* = \\arg \\min_\\lambda f(\\mathbf{x} + \\lambda\\,\\mathbf{v}).\n", "$$\n", "Note that although $f$ is a function of $n$ variables, this minimisation with respect to $\\lambda$ is one dimensional and can be done, for example, using Golden Section Search." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Minimisation in $\\mathbb{R^n}$ : gradient descent\n", "\n", "Gradient descent is one a family of unconstrained optimisation algorithms that can be used when the gradient of the objective function is known or computable.\n", "\n", "* Given $f : \\mathbb{R}^n \\to \\mathbb{R}$ and any point, $\\mathbf{x} \\in \\mathbb{R}^n$, $-\\nabla f(\\mathbf{x})$ points in the direction of steepest decrease of $f$ at $\\mathbf{x}$.\n", "* Idea is to keep going \"downhill\" until a minimum is reached" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Gradient descent algorithm\n", "\n", "Start from a point, $\\mathbf{x}_0$.\n", "\n", "1. Calculate the unit direction vector \n", "$$\\mathbf{v}_n = - \\frac{\\nabla\\,f(\\mathbf{x}_n)}{\\left| \\nabla\\,f(\\mathbf{x}_n)\\right|}.$$\n", "2. Perform line minimisation in direction of $\\mathbf{v}_n$:\n", "$$\n", "\\lambda^* = \\arg \\min_\\lambda f(\\mathbf{x}_n + \\lambda\\,\\mathbf{v}_n)\n", "$$\n", "3. Move to new position $\\mathbf{x}_{n+1} = \\mathbf{x}_n + \\lambda^*\\,\\mathbf{v}_n$. \n", "4. Repeat until $\\left| \\nabla\\,f(\\mathbf{x}_n) \\right| < \\epsilon_\\text{tol}$.\n", "\n", "Normalisation of $\\mathbf{v}_n$ is not strictly necessary but helps keep accuracy as minimum is approached.\n", "\n", "Note: any descent directions will work - gradient is the most efficient." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Gradient descent: example in $\\mathbb{R}^2$\n", "\n", "Consider the function\n", "$$\n", "f(x, y) = \\frac{1}{2}\\left(x^2 + \\gamma y^2 \\right).\n", "$$\n", "The parameter $\\gamma>0$ controls the ellipticity.\n", "The gradient is\n", "$$\n", "\\nabla\\,f(x,y) = (x, \\gamma\\,y)^T.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Gradient descent: example in $\\mathbb{R}^2$ with $\\gamma = 4.0$\n", "\n", "\"array\" " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Gradient descent: example in $\\mathbb{R}^2$ with $\\gamma = 10.0$\n", "\n", "\"array\" " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Gradient descent: Why the zig zags? \n", "\n", "\n", "Remember that $\\mathbf{x}_{n+1}$ is a line minimiser of $f$ from $\\mathbf{x}_n$ in the direction of $\\mathbf{v}_n$. \n", "\n", "Letting $g(\\lambda) = f(\\mathbf{x}_n + \\lambda\\,\\mathbf{v}_n)$, we must have\n", "$$\n", "\\begin{align*}\n", "\\frac{d\\,g}{d\\,\\lambda}(\\lambda^*) &= 0\\\\\n", "\\Rightarrow \\sum_{i=1}^n \\frac{\\partial\\,f}{\\partial\\,x_i}(\\mathbf{x}_n+\\lambda^*\\mathbf{v}_n) \\,v_{n\\,i} &=0\\\\\n", "\\Rightarrow \\nabla\\,f(\\mathbf{x}_{n+1}) \\cdot \\mathbf{v}_n &= 0.\n", "\\end{align*}\n", "$$\n", "\n", "Each step of the gradient descent algorithm is perpendicular to the previous one." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Gradient descent: convergence\n", "\n", "* convergence is exponential but depends on $\\gamma$.\n", "* larger $\\gamma$ is harder - narrower valley requires more zig zags.\n", "* this latter feature means vanilla gradient descent is rarely used in practice.\n", "\n", "\"array\" " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Numerical calculation of derivatives\n", "\n", "Explicit calculation of $\\nabla\\,f(\\mathbf{x})$ is sometimes inconvenient. In such cases it is helpful to be able to calculate derivatives numerically. Several approaches:\n", "\n", "* Finite difference approximation: the workhorse.\n", "* Automatic differentiation: uses dual numbers.\n", "* Spectral methods: use Fourier transforms." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Dual Numbers: basic arithmetic\n", "\n", "Dual numbers are of the form $z = x + \\varepsilon\\, y$ with $x$, $y \\in \\mathbb{R}$ and $\\varepsilon^2 = 0$.\n", "\n", "Addition rule:\n", "$$\n", "\\begin{align*}\n", "z_1 + z_2 & = (x_1 + \\varepsilon\\,y_1) + (x_2 + \\varepsilon\\,y_2)\\\\\n", " & = (x_1 + x_2) + \\varepsilon\\,(y_1+y_2).\n", " \\end{align*}\n", "$$\n", " \n", "Multiplication rule:\n", "$$\n", " \\begin{align*}\n", "z_1 * z_2 & = (x_1 + \\varepsilon\\,y_1) * (x_2 + \\varepsilon\\,y_2) \\\\\n", "& = (x_1 x_2) + \\varepsilon\\,(x_1\\,y_2+x_2\\, y_1).\n", "\\end{align*}\n", "$$\n", "\n", "Dual conjugate:\n", "$$\n", "\\bar{z} = x - \\varepsilon\\, y.\n", "$$\n", "\n", "As with the complex numbers, $z\\,\\bar{z}$ is purely real:\n", "$$\n", "z\\,\\bar{z} = (x + \\varepsilon\\,y) * (x - \\varepsilon\\,y) = x^2. \n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Dual numbers: division and powers\n", "\n", "Division is defined using the conjugate:\n", "$$\n", "\\begin{align*}\n", "\\frac{z_1}{z_2} &= \\frac{z_1\\,\\bar{z}_2}{z_2 \\bar{z}_2}\\\\\n", "&= \\frac{(x_1 + \\varepsilon\\,y_1) (x_2 - \\varepsilon\\,y_2) }{x_2^2}\\\\\n", "&= \\frac{x_1\\,x_2 + \\varepsilon\\,(y_1\\,x_2 - x_1\\,y_2)}{x_2^2}\\\\\n", "&= \\frac{x_1}{x_2} + \\varepsilon\\,\\left( \\frac{y_1\\,x_2 - x_1\\,y_2}{x_2^2}\\right).\n", "\\end{align*}\n", "$$\n", "\n", "Division is not defined for $\\text{Re}(z_2) = 0$ so, unlike $\\mathbb{C}$, the dual numbers do not form a field.\n", "\n", "Powers are obtained using the binomial theorem:\n", "\\begin{align}\n", "\\nonumber (x+ \\varepsilon\\, y)^n & = \\sum_{k=0}^n \\binom{n}{k}\\,x^{n-k}\\, (\\varepsilon\\, y)^k\\\\\n", "\\nonumber &= \\binom{n}{0} x^n + \\binom{n}{1}\\,x^{n-1}\\,\\varepsilon\\, y\\\\\n", "\\label{eq-ADpower} & = x^n + \\varepsilon\\, y\\, n\\,x^{n-1}.\n", "\\end{align}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Dual numbers: automatic differentiation\n", "\n", "Notice that for $f(x) = x^n$,\n", "\n", "$$\n", "f(x+\\varepsilon\\,y) = f(x) + \\varepsilon\\,y\\, f^\\prime(x).\n", "$$\n", "\n", "Evaluating $f(x) = x^n$ at $x + \\varepsilon$ gives the derivative of $f(x)$ as the dual component.\n", "\n", "This is called *automatic differentiation*." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This property trivially extends to polynomial functions and *formally* extends to real analytic functions via their Taylor series:\n", "\n", "$$\n", "\\begin{align*}\n", "&f(x + \\varepsilon\\,y) = \\sum_{k=0}^\\infty \\frac{1}{k!}\\,f^{(k)}(x)\\, \\varepsilon^k\\,y^k\\\\\n", "& = f(x) + \\varepsilon\\,y\\, f^\\prime(x).\n", "\\end{align*}\n", "$$\n", "\n", "Also the chain rule works (check). If $h(x) = f(g(x))$:\n", "\n", "$$\n", "h(x + \\varepsilon) = f(g(x)) + \\varepsilon\\,f^\\prime(g(x))\\,g^\\prime(x).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Dual numbers in Python\n", "\n", "Advantages of automatic differentiation:\n", "* no need to work out complicated analytic formulae for derivatives.\n", "* $\\text{Dual}[z]$ correct to the same precision as $\\text{Real}[z]$.\n", "\n", "Dual arithmetic is supported in Python, with [num_dual](https://pypi.org/project/num-dual/) discussed below." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2 + [7]ε" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Install e.g. via \"python3 -m pip install num_dual\" (or pip install)\n", "\n", "from num_dual import Dual64\n", "\n", "z1 = Dual64(2.0,3.0)\n", "z2 = Dual64(1.0,2.0)\n", "z1*z2" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The value of the function at the target point is 4.0\n", "Its first derivative is 8.0\n", "Its second derivative is 27.0\n", "Its third derivative is 18.0\n" ] } ], "source": [ "# f'(x) = x^3, f''(x) = 3*x^2, f'''(x) = 6*x\n", "def f(x):\n", " return 0.25*x**4\n", "\n", "# Evaluate at 2.0 using the dual part set to 1.0\n", "x = Dual64(2.0, 1.0)\n", "\n", "result = f(x)\n", "\n", "print(\"The value of the function at the target point is\", result.value)\n", "print(\"Its first derivative is\", result.first_derivative)\n", "\n", "# If we want to go further ...\n", "from num_dual import derive3\n", "\n", "x = derive3(3.0)\n", "result = f(x)\n", "\n", "print(\"Its second derivative is\", result.second_derivative)\n", "print(\"Its third derivative is\", result.third_derivative)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Automatic differentiation for multivariate functions\n", "\n", "In the formulation presented above, automatic differentiation of multivariate functions requires a separate function evaluation for each partial derivative.\n", "\n", "For a multivariate function $f : \\mathbb{R}^n \\to \\mathbb{R}$, to calculate the partial derivative with respect to $x_i$, we have to add $\\varepsilon$ to the $i^\\text{th}$ component of the argument and then take the dual part.\n", "For example, in $\\mathbb{R}^2$ we have\n", "$$\n", "\\begin{align*}\n", "f(x + \\varepsilon, y) &= f(x,y) + \\varepsilon\\,\\frac{\\partial f}{\\partial x} (x,y)\\\\\n", "f(x, y + \\varepsilon) &= f(x,y) + \\varepsilon\\,\\frac{\\partial f}{\\partial y} (x,y).\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Minimisation of sums\n", "\n", "Problems in statistical inference, ML and data science frequently produce optimisation problems that involve minimising a sum:\n", "$$\n", "\\min_\\beta f(\\beta) = \\min_\\beta \\frac{1}{n}\\sum_{i=1}^n L(\\beta, \\mathbf{x}_i).\n", "$$\n", "Often $\\beta$ represents some model parameters, $\\mathbf{x}_i$ represents the $i^\\text{th}$ observation in the training data set and $L$ is a loss function of some kind.\n", "\n", "Can be solved with GD by calculating the gradient with respect to $\\beta$:\n", "$$\n", "(\\nabla f(\\beta))_j = \\frac{1}{n}\\sum_{i=1}^n \\frac{\\partial\\, L}{\\partial\\, \\beta_j}(\\beta, \\mathbf{x}_i)\n", "$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Stochastic Gradient Descent\n", "\n", "If the training data is large ($n \\gg 1$), then calculating the gradient can become very expensive. Stochastic Gradient Descent addresses this issue.\n", "\n", "Idea is to estimate the gradients as:\n", "\n", "$$\n", "(\\nabla f(\\beta))_j \\approx (\\tilde{\\nabla} f(\\beta))_j = \\frac{\\partial\\, L}{\\partial\\, \\beta_j}(\\beta, \\mathbf{\\tilde{x}}_{i})\n", "$$\n", "\n", "where $\\tilde{\\mathbf{x}}_i$ is the $i^\\text{th}$ training data point from the training data set having first undergone a random shuffle.\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Stochastic Gradient Descent : learning rate\n", "\n", "With SGD, we typically do not perform full line minimisations (due to sampling noise). Instead the basic update rule is\n", "$$\n", "\\mathbf{\\beta}_{n+1} = \\mathbf{\\beta}_n - \\zeta_n\\,(\\tilde{\\nabla} f(\\beta_n)). \n", "$$\n", "$\\zeta_n$ is a (decreasing) function of $n$ often called the *learning rate*. For example,\n", "\n", "$$\n", "\\zeta_n = \\frac{\\zeta_0}{1+\\zeta_1\\,n},\n", "$$\n", "\n", "where $\\zeta_0$ and $\\zeta_1$ are *hyperparameters*. Note that $\\sum_{n=1}^\\infty \\zeta_n$ is a divergent series. Prevents SGD from stalling.\n", "\n", "Due to sampling noise, SGD only converges to a \"noise ball\" around the minimum of $f(\\beta)$ rather than to the true minimum." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Logistic regression\n", "\n", "Binary variable, $Y \\in \\left\\{0, 1\\right\\}$, explanatory variable, $X \\in \\mathbb{R}$. \n", "\n", "Example: how is eye strain ($Y$) related to the number of hours per week ($X$) spent on MS Teams?\n", "\n", "\"array\" " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Assume $Y = \\left\\{Y_i, i=1\\ldots n \\right\\}$ are id Bernoulli with parameters, $p_i$, depending on $X = \\left\\{X_i, i=1\\ldots n \\right\\}$:\n", "$$\n", "\\mathbb{P}(Y | X) = \\prod_{i=1}^n p(X_i)^{Y_i} \\, (1 -p(X_i))^{1-Y_i}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Logistic regression\n", "\n", "How should the Bernoulli parameter, $p$, vary with $X$? \n", "Assume *log odds* are a linear function of $X$:\n", "$$\n", "\\log\\left(\\frac{p(X)}{1-p(X)} \\right) = \\beta_0 + \\beta_1\\,X = \\beta \\cdot \\mathbf{X}\n", "$$\n", "where $\\mathbf{X} = \\left(1, X\\right)$ and $\\beta = \\left(\\beta_0, \\beta_1\\right)$.\n", "\n", "Solve for $p(X)$:\n", "$$\n", "p(X) = \\sigma(\\beta\\cdot\\mathbf{X}),\n", "$$\n", "where $\\sigma(x)$ is the sigmoid function:\n", "$$\n", "\\sigma(x) = \\frac{1}{1+\\exp(-x)}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Logistic regression as an optimisation problem\n", "\n", "Obviously cannot use least squares to select optimal values of $\\beta$. \n", "\n", "Instead maximise, $\\mathbb{P}(Y | X) $ - *likelihood of $Y$ given $X$* - with respect to $\\beta$. \n", "\n", "In fact easier to minimise $-\\log \\mathbb{P}(Y | X)$ because it turns task into a sum minimisation: \n", "\n", "$$\n", "\\beta_* = \\arg \\min_{\\beta} \\sum_{i=1}^n -Y_i\\,\\log(\\sigma(\\beta\\cdot\\mathbf{X}_i)) - (1 -Y_i)\\,\\log(1-\\sigma(\\beta\\cdot\\mathbf{X}_i)).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Logistic regression via gradient descent\n", "\n", "We can solve this sum minimisation problem with GD or SGD. First we need the gradient of the objective function,\n", "$$\n", "L(\\beta) = \\sum_{i=1}^n -Y_i\\,\\log(\\sigma(\\beta\\cdot\\mathbf{X}_i)) - (1 -Y_i)\\,\\log(1-\\sigma(\\beta\\cdot\\mathbf{X}_i)).\n", "$$\n", "Direct calculation gives (check)\n", "$$\n", "(\\nabla L (\\beta))_k = \\sum_{i=1}^n (\\sigma(\\beta\\cdot\\mathbf{X}_i)) - Y_i)\\,(\\mathbf{X}_i)_k\n", "$$\n", "where $k\\in\\left\\{0,1\\right\\}$ and $(\\mathbf{X}_i)_k$ is the k$^\\text{th}$ component of $\\mathbf{X}_i$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implementation: \n", "\n", "1. Generate some test data\n", "2. Define the objective function\n", "\n", "Need negative log likelihood of the data $(\\mathbf{X}, \\mathbf{Y})$ for a given value of $\\beta$. \n", "First calculate the components of the sum as separate lists and then assemble them at end.\n", "\n", "3. Define the gradient of the objective function\n", "4. Define the partial gradient (for SGD)\n", "5. Run gradient descent\n", "6. Run stochastic gradient descent" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot of L(β)\n", "\n", "\"array\" " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Comparison of results\n", "\n", "\"array\" " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Constrained nonlinear optimisation\n", "\n", "Given $f(\\mathbf{x}) : \\mathbb{R}^n \\to \\mathbb{R}$, we can combine what we have learned about gradient descent with what we will (soon) learn about linear programming to solve constrained nonlinear optimisation problems of the form:\n", "$$\n", "\\min_{\\mathbf{x}\\in C} f(\\mathbf{x}), \n", "$$\n", "where $C \\subset \\mathbb{R}^n$ is defined by a set of linear inequalities.\n", "\n", "Clearly GD alone is insufficient since $\\mathbf{x}_n$ could leave the feasible set, $C$.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Constrained nonlinear optimisation: Frank-Wolfe algorithm\n", "\n", "Idea is to linearise $f(\\mathbf{x})$ about $\\mathbf{x}_n$ at each iteration of GD and solve the associated LP to obtain an $\\mathbf{x}_{n+1}$, which cannot leave $C$ by construction.\n", "\n", "Given current point, $\\mathbf{x}_n$, \n", "\n", "$$\n", "f(\\mathbf{x}) = f(\\mathbf{x}_n)+(\\mathbf{x}-\\mathbf{x}_n) \\cdot \\nabla\\,f(\\mathbf{x}_n) + \\mathcal{O}(\\left| \\mathbf{x}-\\mathbf{x}_n \\right|^2).\n", "$$\n", "\n", "Neglecting the $\\mathcal{O}(\\left| \\mathbf{x}-\\mathbf{x}_n \\right|^2)$ terms, minimisation of $f(\\mathbf{x})$, is now a LP:\n", "\n", "$$\n", "x^* = \\arg \\min_{\\mathbf{x}\\in C}\\ \\left[ \\mathbf{x} \\cdot \\nabla\\,f(\\mathbf{x}_n)\\right]\n", "$$\n", "\n", "Start from a point, $\\mathbf{x}_0 \\in C$. At each iteration, $n$:\n", "\n", "1. Solve the linearised problem\n", "$$\n", "x^* = \\arg \\min_{\\mathbf{x}\\in C}\\ \\left[\\mathbf{x} \\cdot \\nabla\\,f(\\mathbf{x}_n)\\right].\n", "$$\n", "2. Perform *bounded* line minimisation of $f(\\mathbf{x})$ from $\\mathbf{x}_n$ in the direction of $\\mathbf{x}^*-\\mathbf{x}_n$: \n", "$$\n", "\\lambda^* = \\arg \\min_{\\lambda \\in \\left[0,1\\right]} f \\left( \\mathbf{x}_n + \\lambda\\,(\\mathbf{x}^*-\\mathbf{x}_n)\\right).\n", "$$\n", "3. Update\n", "$$\n", "\\mathbf{x}_{n+1} = \\mathbf{x}_n + \\lambda^* \\,(\\mathbf{x}^*-\\mathbf{x}_n).\n", "$$\n", "4. Repeat until desired tolerance is reached." ] } ], "metadata": { "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Python 3", "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.8.5" } }, "nbformat": 4, "nbformat_minor": 2 }