From 655f5adc7254e5e3286a95ddeb7d0fa0c6110b15 Mon Sep 17 00:00:00 2001 From: Brian Delhaisse Date: Tue, 17 Dec 2019 00:09:34 +0100 Subject: [PATCH] add few tutorials (under construction) --- tutorials/README.md | 6 + tutorials/machine_learning/GP.ipynb | 808 +++++++++++++ tutorials/machine_learning/README.md | 25 + tutorials/math/Linear-Algebra.ipynb | 487 ++++++++ tutorials/math/PCA.ipynb | 1144 +++++++++++++++++++ tutorials/math/README.md | 14 + tutorials/python/README.md | 65 ++ tutorials/robotics/README.md | 27 + tutorials/robotics/cartpole.ipynb | 39 + tutorials/robotics/dynamics.ipynb | 39 + tutorials/robotics/kinematics.ipynb | 39 + tutorials/robotics/mass-spring-damper.ipynb | 215 ++++ 12 files changed, 2908 insertions(+) create mode 100644 tutorials/README.md create mode 100644 tutorials/machine_learning/GP.ipynb create mode 100644 tutorials/machine_learning/README.md create mode 100644 tutorials/math/Linear-Algebra.ipynb create mode 100644 tutorials/math/PCA.ipynb create mode 100644 tutorials/math/README.md create mode 100644 tutorials/python/README.md create mode 100644 tutorials/robotics/README.md create mode 100644 tutorials/robotics/cartpole.ipynb create mode 100644 tutorials/robotics/dynamics.ipynb create mode 100644 tutorials/robotics/kinematics.ipynb create mode 100644 tutorials/robotics/mass-spring-damper.ipynb diff --git a/tutorials/README.md b/tutorials/README.md new file mode 100644 index 0000000..5f73a06 --- /dev/null +++ b/tutorials/README.md @@ -0,0 +1,6 @@ +## Tutorials + +UNDER CONSTRUCTION + +This folder contains tutorials about robotics and machine learning in the form of python notebooks and/or webpages. + diff --git a/tutorials/machine_learning/GP.ipynb b/tutorials/machine_learning/GP.ipynb new file mode 100644 index 0000000..10263ca --- /dev/null +++ b/tutorials/machine_learning/GP.ipynb @@ -0,0 +1,808 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Gaussian Process\n", + "\n", + "In this tutorial, we expose what gaussian processes are, and how to use the [GPy library](http://sheffieldml.github.io/GPy/). We first provide a gentle reminder about Gaussian distributions and their properties." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Import all the important libraries\n", + "import math\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from scipy.stats import norm\n", + "from scipy.stats import multivariate_normal\n", + "from mpl_toolkits.mplot3d import Axes3D\n", + "import matplotlib.mlab as mlab\n", + "from ipywidgets import widgets as wg\n", + "from matplotlib import cm\n", + "\n", + "import GPy\n", + "#%matplotlib inline\n", + "%matplotlib notebook" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1D Gaussian distribution" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot\n", + "def plot_gaussian(mu=0, sigma=1):\n", + " x = np.linspace(-3, 3, 100)\n", + " plt.plot(x, norm.pdf(x, mu, sigma))\n", + " plt.xlabel('x')\n", + " plt.ylabel('p(x)')\n", + "\n", + "wg.interact(plot_gaussian, mu=(-2,2,0.1), sigma=(-2,2,0.1))\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Multivariate Gaussian distribution (2D)\n", + "\n", + "The multivariable Gaussian distribution is a generalization of the Gaussian distribution to vectors. See [wikipedia](https://en.wikipedia.org/wiki/Multivariate_normal_distribution) for more info." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# moments\n", + "mu = np.array([0,0])\n", + "Sigma = np.array([[1,0], \n", + " [0,1]])\n", + "Sigma1 = np.array([[1,0.5],\n", + " [0.5,1]])\n", + "Sigma2 = np.array([[1,-0.5],\n", + " [-0.5,1]])\n", + "Sigmas = [Sigma, Sigma1, Sigma2]\n", + "\n", + "pts = []\n", + "for S in Sigmas:\n", + " pts.append(np.random.multivariate_normal(mu, S, 1000).T)\n", + "\n", + "# Plotting\n", + "width = 16\n", + "height = 4\n", + "plt.figure(figsize=(width, height))\n", + "\n", + "# make plot\n", + "for i in range(len(Sigmas)):\n", + " plt.subplot(1,3,i+1)\n", + " plt.title('Plot '+str(i+1))\n", + " plt.ylim(-4,4)\n", + " plt.xlim(-4,4)\n", + " plt.xlabel('x1')\n", + " plt.ylabel('x2')\n", + " plt.plot(pts[i][0], pts[i][1],'o')\n", + " #plt.scatter(pts[i][0], pts[i][1])\n", + " \n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The 1st plot above is described by:\n", + "\\begin{equation}\n", + " \\left[ \\begin{array}{c} x_1\\\\x_2 \\end{array}\\right] \\sim \\mathcal{N} \\left(\\left[ \\begin{array}{c} 0\\\\0 \\end{array}\\right], \\left[ \\begin{array}{cc} 1 & 0\\\\ 0 & 1 \\end{array}\\right] \\right)\n", + "\\end{equation}\n", + "\n", + "The 2nd plot is given by:\n", + "\\begin{equation}\n", + " \\left[ \\begin{array}{c} x_1\\\\x_2 \\end{array}\\right] \\sim \\mathcal{N}\\left(\\left[ \\begin{array}{c} 0\\\\0 \\end{array}\\right], \\left[ \\begin{array}{cc} 1 & 0.5\\\\ 0.5 & 1 \\end{array}\\right]\\right)\n", + "\\end{equation}\n", + "\n", + "Finally, the 3rd plot is given by:\n", + "\\begin{equation}\n", + " \\left[ \\begin{array}{c} x_1\\\\x_2 \\end{array}\\right] \\sim \\mathcal{N}\\left(\\left[ \\begin{array}{c} 0\\\\0 \\end{array}\\right], \\left[ \\begin{array}{cc} 1 & -0.5\\\\ -0.5 & 1 \\end{array}\\right]\\right)\n", + "\\end{equation}\n", + "\n", + "The covariance (and the dot product) measures the similarity.\n", + "\n", + "For the 2nd and 3rd plots, $x_1$ is **correlated** with $x_2$, i.e. knowing $x_1$ gives us information about $x_2$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Joint distribution $p(x_1,x_2)$\n", + "\n", + "The joint distribution $p(x_1, x_2)$ is given by:\n", + "\n", + "\\begin{equation}\n", + " \\left[ \\begin{array}{c} x_1\\\\x_2 \\end{array}\\right] \\sim \\mathcal{N}\\left(\\left[ \\begin{array}{c} \\mu_1 \\\\ \\mu_2 \\end{array}\\right], \\left[ \\begin{array}{cc} \\Sigma_{11} & \\Sigma_{12} \\\\ \\Sigma_{21} & \\Sigma_{22} \\end{array}\\right] \\right) = \\mathcal{N}(\\pmb{\\mu}, \\pmb{\\Sigma})\n", + "\\end{equation}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Reference: http://stackoverflow.com/questions/38698277/plot-normal-distribution-in-3d\n", + "\n", + "# moments\n", + "mu = np.array([0,0])\n", + "Sigma = np.array([[1,0], [0,1]])\n", + "\n", + "# Create grid and multivariate normal\n", + "step = 500\n", + "bound = 10\n", + "x = np.linspace(-bound,bound,step)\n", + "y = np.linspace(-bound,bound,step)\n", + "X, Y = np.meshgrid(x,y)\n", + "pos = np.empty(X.shape + (2,))\n", + "pos[:, :, 0] = X; pos[:, :, 1] = Y\n", + "pdf = multivariate_normal(mu, Sigma).pdf(pos)\n", + "\n", + "# Plot\n", + "fig = plt.figure(figsize=plt.figaspect(0.5)) # Twice as wide as it is tall.\n", + "\n", + "# 1st subplot (3D)\n", + "ax = fig.add_subplot(1, 2, 1, projection='3d')\n", + "ax.plot_surface(X, Y, pdf, cmap='viridis', linewidth=0)\n", + "ax.set_xlabel('x1')\n", + "ax.set_ylabel('x2')\n", + "ax.set_zlabel('p(x1, x2)')\n", + "\n", + "# 2nd subplot (2D)\n", + "ax = fig.add_subplot(1, 2, 2)\n", + "ax.contourf(x, y, pdf)\n", + "#ax.colorbar()\n", + "ax.set_xlabel('x1')\n", + "ax.set_ylabel('x2')\n", + "\n", + "fig.tight_layout()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Normalization\n", + "\n", + "In order to be a valid probability distribution, the volume under the surface should equal to 1.\n", + "\n", + "\\begin{equation}\n", + " \\int \\int p(x_1,x_2) dx_1 dx_2 = 1\n", + "\\end{equation}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# p(x1,x2) = pdf\n", + "# dx = 2.*bound/step\n", + "# dx1 dx2 = (2.*bound/step)**2\n", + "print(\"Summation: {}\".format((2.*bound/step)**2 * pdf.sum()))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Conditional distribution $p(x_2|x_1)$\n", + "\n", + "What is the mean $\\mu_{2|1}$ and the variance $\\Sigma_{2|1}$ of the conditional distribution $p(x_2|x_1) = \\mathcal{N}(\\mu_{2|1}, \\Sigma_{2|1})$?\n", + "\n", + "We know the mean $\\pmb{\\mu}$ and the covariance $\\pmb{\\Sigma}$ of the joint distribution $p(x_1,x_2)$. Using the [Schur complement](https://en.wikipedia.org/wiki/Schur_complement), we obtain:\n", + "\n", + "\\begin{align}\n", + " \\mu_{2|1} &= \\mu_{2} + \\Sigma_{21}\\Sigma_{22}^{-1}(x_2 - \\mu_2) \\\\\n", + " \\Sigma_{2|1} &= \\Sigma_{22} - \\Sigma_{21}\\Sigma_{22}^{-1}\\Sigma_{12}\n", + "\\end{align}\n", + "\n", + "For the demo, check Murphy's book \"Machine Learning: A Probabilistic Perspective\", section 4.3.4" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=plt.figaspect(0.5)) # Twice as wide as it is tall.\n", + "\n", + "x1_value = 0\n", + "z_max = pdf.max()\n", + "\n", + "# 1st subplot\n", + "ax = fig.add_subplot(1, 2, 1, projection='3d')\n", + "ax.plot_surface(X, Y, pdf, cmap='viridis', linewidth=0)\n", + "ax.set_xlabel('x1')\n", + "ax.set_ylabel('x2')\n", + "ax.set_zlabel('p(x1, x2)')\n", + "y1 = np.linspace(-bound,bound,2)\n", + "z = np.linspace(0,z_max,2)\n", + "Y1, Z = np.meshgrid(y1,z)\n", + "ax.plot_surface(x1_value, Y1, Z, color='red', alpha=0.2)\n", + "#cset = ax.contourf(X, Y, pdf, zdir='x', offset=-bound, cmap=cm.coolwarm)\n", + "\n", + "# 2nd subplot\n", + "ax = fig.add_subplot(1, 2, 2)\n", + "ax.plot(x, pdf[step//2 + x1_value*step//(2*bound)])\n", + "ax.set_xlabel('x2')\n", + "ax.set_ylabel('p(x2|x1)')\n", + "\n", + "fig.tight_layout()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Marginal distribution $p(x_1)$ and $p(x_2)$\n", + "\n", + "\\begin{align}\n", + " p(x_1) &= \\int p(x_1, x_2) dx_2 = \\mathcal{N}(\\mu_1, \\Sigma_{11}) \\\\\n", + " p(x_2) &= \\int p(x_1, x_2) dx_1 = \\mathcal{N}(\\mu_2, \\Sigma_{22})\n", + "\\end{align}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=plt.figaspect(0.5))\n", + "plt.subplot(1,2,1)\n", + "plt.title('By summing')\n", + "dx = 2. * bound / step\n", + "plt.plot(x, pdf.sum(0) * dx, color='blue')\n", + "plt.xlabel('x2')\n", + "plt.ylabel('p(x2)')\n", + "\n", + "plt.subplot(1,2,2)\n", + "plt.title('by using the normal distribution')\n", + "plt.plot(x, norm.pdf(x, mu[1], Sigma[1,1]), color='red')\n", + "plt.xlabel('x2')\n", + "plt.ylabel('p(x2)')\n", + "\n", + "fig.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Gaussian Processes (GPs)\n", + "\n", + "A Gaussian process is a Gaussian distribution over functions. That is, it is a generalization of the multivariable Gaussian distribution to infinite vectors.\n", + "\n", + "It will become clearer with an example." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x = [0.5,0.8,1.4]\n", + "f = [1,2,6]\n", + "\n", + "plt.plot(x,f,'o')\n", + "for i in range(len(x)):\n", + " plt.annotate('f'+str(i+1), (x[i],f[i]))\n", + "plt.xlim(0,2)\n", + "plt.ylim(0,6.5)\n", + "plt.ylabel('f(x)')\n", + "plt.xlabel('x')\n", + "plt.xticks(x, ['x'+str(i+1) for i in range(len(x))])\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\\begin{align}\n", + " \\left[ \\begin{array}{c} f_1 \\\\ f_2 \\\\ f_3 \\end{array}\\right] \n", + " &\\sim \\mathcal{N}\\left( \\left[ \\begin{array}{c} 0 \\\\ 0 \\\\ 0 \\end{array}\\right], \\left[ \\begin{array}{ccc} K_{11} & K_{12} & K_{13} \\\\ K_{21} & K_{22} & K_{23} \\\\ K_{31} & K_{32} & K_{33} \\end{array}\\right] \\right) \\\\\n", + " &\\sim \\mathcal{N}\\left( \\left[ \\begin{array}{c} 0 \\\\ 0 \\\\ 0 \\end{array}\\right], \\left[ \\begin{array}{ccc} 1 & 0.7 & 0.2 \\\\ 0.7 & 1 & 0.6 \\\\ 0.2 & 0.6 & 1 \\end{array}\\right] \\right)\n", + "\\end{align}\n", + "\n", + "Similarity measure: $K_{ij} = \\exp(- ||x_i - x_j||^2) = \\left\\{ \\begin{array}{ll} 0 & ||x_i - x_j|| \\rightarrow \\infty \\\\ 1 & x_i = x_j \\end{array} \\right.$\n", + "\n", + "Prediction (noiseless GP regression): given data $\\mathcal{D} = \\{(x_1,f_1), (x_2,f_2), (x_3,f_3)\\}$, and new point $x_*$ (e.g. $x_*$=1.4), what is the value of $f_*$?\n", + "\n", + "\\begin{equation}\n", + " \\pmb{f} \\sim \\mathcal{N}(\\pmb{0}, \\pmb{K}) \\qquad \\mbox{and} \\qquad f_* \\sim \\mathcal{N}(0, K(x_*,x_*)) = \\mathcal{N}(0, K_{**})\n", + "\\end{equation}\n", + "\n", + "In this case, $K_{**} = K(x_*,x_*) = \\exp(- ||x_* - x_*||^2) = 1$.\n", + "\n", + "Now, we can write:\n", + "\n", + "\\begin{equation}\n", + " \\left[ \\begin{array}{c} \\pmb{f} \\\\ f_* \\end{array} \\right] \\sim \\mathcal{N}\\left(\\pmb{0}, \\left[\\begin{array}{cc} \\left[ \\begin{array}{ccc} K_{11} & K_{12} & K_{13} \\\\ K_{21} & K_{22} & K_{23} \\\\ K_{31} & K_{32} & K_{33} \\end{array} \\right] & \\left[ \\begin{array}{c} K_{1*} \\\\ K_{2*} \\\\ K_{3*} \\end{array} \\right] \\\\ \\left[ \\begin{array}{ccc} K_{*1} & K_{*2} & K_{*3} \\end{array} \\right] & \\left[\\begin{array}{c} K_{**} \\end{array} \\right] \\end{array} \\right] \\right) = \\mathcal{N}\\left(\\pmb{0}, \\left[\\begin{array}{cc} \\pmb{K} & \\pmb{K}_* \\\\ \\pmb{K}_*^T & \\pmb{K}_{**} \\end{array}\\right]\\right)\n", + "\\end{equation}\n", + "\n", + "Using the formula for the conditional probability $p(f_*|f)$, we have:\n", + "\n", + "\\begin{align}\n", + " \\mu_* &= \\mathbb{E}[f_*] = \\pmb{K}_*^T \\pmb{K}^{-1}\\pmb{f} \\\\\n", + " c_* &= K_{**} - \\pmb{K}_*^T \\pmb{K}^{-1}\\pmb{K}_*\n", + "\\end{align}\n", + "\n", + "We can thus predict the mean $\\mu_*$ and the variance $c_*$ for the test point $x_*$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x = [0.5,0.8,1.4]\n", + "f = [1,2,6]\n", + "\n", + "x_new = 1.3\n", + "f_new = 5.2\n", + "\n", + "plt.plot(x+[x_new],f+[f_new],'o')\n", + "for i in range(len(x)):\n", + " plt.annotate('f'+str(i+1), (x[i],f[i]))\n", + "plt.errorbar(x_new, f_new, yerr=1)\n", + "plt.annotate('f*', (x_new+0.02, f_new))\n", + "plt.xlim(0,2)\n", + "plt.ylim(0,6.5)\n", + "plt.ylabel('f(x)')\n", + "plt.xlabel('x')\n", + "plt.xticks(x+[x_new], ['x'+str(i+1) for i in range(len(x))]+['x*'])\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Generalization\n", + "\n", + "A GP defines a distribution over functions $p(f)$ (i.e. it is the joint distribution over all the infinite function values).\n", + "\n", + "Definition: $p(f)$ is a GP if for any finite subset $\\{x_1,...,x_n\\} ⊂ X$, the marginal distribution over that finite subset $p(f)$ has a multivariate Gaussian distribution.\n", + "\n", + "Prior on $f$:\n", + "\\begin{equation}\n", + " \\pmb{f}|\\pmb{x} \\sim \\mathcal{GP}(\\pmb{\\mu}(\\pmb{x}), \\pmb{K}(\\pmb{x}, \\pmb{x}))\n", + "\\end{equation}\n", + "with\n", + "\\begin{align*}\n", + " \\pmb{\\mu}(\\pmb{x}) &= \\mathbb{E}_f \\lbrack \\pmb{x} \\rbrack \\\\\n", + " k(\\pmb{x}, \\pmb{x'}) &= \\mathbb{E}_f \\lbrack (\\pmb{x} - \\pmb{\\mu}(\\pmb{x})) (\\pmb{x'} - \\pmb{\\mu}(\\pmb{x'})) \\rbrack\n", + "\\end{align*}\n", + "\n", + "Often written as:\n", + "\\begin{equation}\n", + " \\pmb{f} \\sim \\mathcal{GP}(\\pmb{0}, \\pmb{K})\n", + "\\end{equation}\n", + "\n", + "Concretely, assume $\\pmb{x} \\in \\mathbb{R}^{50}$, then $\\pmb{K}(\\pmb{x}, \\pmb{x}) \\in \\mathbb{R}^{50 \\times 50}$, then $\\pmb{f} \\sim \\mathcal{GP}(\\pmb{0}, \\pmb{K})$ means:\n", + "\n", + "\\begin{equation}\n", + " \\left[ \\begin{array}{c} f_1 \\\\ \\vdots \\\\ f_{50} \\end{array}\\right] := \\left[ \\begin{array}{c} f(x_1) \\\\ \\vdots \\\\ f(x_{50}) \\end{array}\\right] \\sim \\mathcal{N}\\left( \\left[ \\begin{array}{c} 0 \\\\ \\vdots \\\\ 0 \\end{array}\\right], \\left[ \\begin{array}{ccc} k(x_1,x_1) & \\cdots & k(x_1, x_{50}) \\\\ \\vdots & \\ddots & \\vdots \\\\ k(x_{50},x_1) & \\cdots & k(x_{50}, x_{50}) \\end{array} \\right] \\right)\n", + "\\end{equation}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### RBF kernel\n", + "\n", + "Let's choose a RBF (a.k.a Squared Exponential, Gaussian) kernel:\n", + "\n", + "\\begin{equation}\n", + "\\pmb{K} = \\left[ \\begin{array}{ccc} k(x_1,x_1) & \\cdots & k(x_1, x_d) \\\\ \\vdots & \\ddots & \\vdots \\\\ k(x_d,x_1) & \\cdots & k(x_d, x_d) \\end{array} \\right]\n", + "\\end{equation}\n", + "with\n", + "\\begin{equation}\n", + " k(x_i, x_j) = \\alpha^2 \\exp \\left( - \\frac{(x_i - x_j)^2}{2l} \\right) \\qquad \\mbox{ and hyperparameters } \\pmb{\\Phi} = \\left\\{ \\begin{array}{l} \\alpha \\mbox{: amplitude} \\\\ l \\mbox{: the lengthscale} \\end{array} \\right.\n", + "\\end{equation}\n", + "\n", + "This function $k$ is infinitely differentiable." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Reference: https://www.youtube.com/watch?v=4vGiHC35j9s&t=51s\n", + "\n", + "# Hyperparameters\n", + "alpha = 1\n", + "l = 2\n", + "\n", + "# Parameters\n", + "n = 50 # nb of points\n", + "n_func = 10 # nb of fct to draw\n", + "x_bound = 5 # bound on the x axis\n", + "\n", + "def RBF_kernel(a,b):\n", + " sqdist = np.sum(a**2,1).reshape(-1,1) + np.sum(b**2,1) - 2*np.dot(a,b.T)\n", + " return alpha**2 * np.exp(-1/l * sqdist)\n", + "\n", + "n = 50\n", + "X = np.linspace(-x_bound, x_bound, n).reshape(-1,1)\n", + "K = RBF_kernel(X, X) # dim(K) = n x n\n", + "\n", + "L = np.linalg.cholesky(K + 1e-6 * np.eye(n))\n", + "f_prior = np.dot(L, np.random.normal(size=(n, n_func)))\n", + "\n", + "# Plotting\n", + "width = 16\n", + "height = 4\n", + "plt.figure(figsize=(width, height))\n", + "\n", + "# plot f_prior\n", + "plt.subplot(1,3,1)\n", + "plt.title('GP: prior on f')\n", + "plt.plot(X, f_prior)\n", + "plt.plot(X, f_prior.mean(1), linewidth=3, color='black')\n", + "plt.ylabel('f(x)')\n", + "plt.xlabel('x')\n", + "\n", + "# plot Kernel\n", + "plt.subplot(1,3,2)\n", + "plt.title('Kernel matrix')\n", + "plt.pcolor(K[::-1])\n", + "plt.colorbar()\n", + "\n", + "plt.subplot(1,3,3)\n", + "plt.title('Kernel function')\n", + "plt.plot(X, RBF_kernel(X, np.array([[1.0]])))\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Kernel (prior knowledge)\n", + "\n", + "By choosing a specific kernel, we can incorporate prior knowledge that we have about the function $f$, such as, if the function is:\n", + "* periodic\n", + "* smooth\n", + "* symmetric\n", + "* etc.\n", + "\n", + "The hyperparameters for each kernel are also very intuitive/interpretable.\n", + "\n", + "Note: kernels can be combined!\n", + "\n", + "Indeed, if $k(x,y)$, $k_1(x,y)$ and $k_2(x,y)$ are valid kernels then:\n", + "* $\\alpha k(x,y) $ with $\\alpha \\geq 0$\n", + "* $k_1(x,y) + k_2(x,y)$\n", + "* $k_1(x,y) k_2(x,y)$\n", + "* $p(k(x,y))$ with $p$ being a polynomial function with non-negative coefficients\n", + "* $exp(k(x,y))$\n", + "* $f(x) k(x,y) \\overline{f(y)}$ with $\\overline{f} = $ complex conjugate\n", + "* $k(\\phi(x),\\phi(y))$\n", + "\n", + "are all valid kernels!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Periodic Exponential kernel" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "variance = 1.\n", + "lengthscale = 1.\n", + "period = 2.*np.pi\n", + "\n", + "#K = periodic_kernel(X, X) # dim(K) = n x n\n", + "kern = GPy.kern.PeriodicExponential(variance=variance, lengthscale=lengthscale, period=period)\n", + "K1 = kern.K(X)\n", + "\n", + "L = np.linalg.cholesky(K1 + 1e-6 * np.eye(n))\n", + "f_prior = np.dot(L, np.random.normal(size=(n, 1)))\n", + "\n", + "# Plotting\n", + "width = 16\n", + "height = 4\n", + "plt.figure(figsize=(width, height))\n", + "\n", + "# plot f_prior\n", + "plt.subplot(1,3,1)\n", + "plt.title('GP: prior on f')\n", + "plt.plot(X, f_prior)\n", + "plt.plot(X, f_prior.mean(1), linewidth=3, color='black')\n", + "plt.ylabel('f(x)')\n", + "plt.xlabel('x')\n", + "\n", + "# plot Kernel\n", + "plt.subplot(1,3,2)\n", + "plt.title('Kernel matrix')\n", + "plt.pcolor(K1[::-1])\n", + "plt.colorbar()\n", + "\n", + "plt.subplot(1,3,3)\n", + "plt.title('Kernel function')\n", + "plt.plot(X, kern.K(X, np.array([[1.0]])))\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### addition and multiplication of 2 kernels (SE and PE)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "K_add = K + K1\n", + "\n", + "L = np.linalg.cholesky(K_add + 1e-6 * np.eye(n))\n", + "f_prior = np.dot(L, np.random.normal(size=(n, n_func)))\n", + "\n", + "# Plotting\n", + "width = 16\n", + "height = 8\n", + "plt.figure(figsize=(width, height))\n", + "\n", + "# plot f_prior\n", + "plt.subplot(2,2,1)\n", + "plt.title('GP: prior on f with K_add')\n", + "plt.plot(X, f_prior)\n", + "plt.plot(X, f_prior.mean(1), linewidth=3, color='black')\n", + "plt.ylabel('f(x)')\n", + "plt.xlabel('x')\n", + "\n", + "# plot Kernel\n", + "plt.subplot(2,2,2)\n", + "plt.title('Kernel matrix: K_add')\n", + "plt.pcolor(K_add[::-1])\n", + "plt.colorbar()\n", + "\n", + "K_prod = K * K1\n", + "\n", + "L = np.linalg.cholesky(K_prod + 1e-6 * np.eye(n))\n", + "f_prior = np.dot(L, np.random.normal(size=(n, n_func)))\n", + "\n", + "# plot f_prior\n", + "plt.subplot(2,2,3)\n", + "plt.title('GP: prior on f with K_prod')\n", + "plt.plot(X, f_prior)\n", + "plt.plot(X, f_prior.mean(1), linewidth=3, color='black')\n", + "plt.ylabel('f(x)')\n", + "plt.xlabel('x')\n", + "\n", + "# plot Kernel\n", + "plt.subplot(2,2,4)\n", + "plt.title('Kernel matrix: K_prod')\n", + "plt.pcolor(K_prod[::-1])\n", + "plt.colorbar()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### GP Posterior\n", + "\n", + "Given $\\mathcal{D}=\\{(x_i, y_i)\\}_{i=1}^{i=N} = (\\pmb{X}, \\pmb{y})$, we have:\n", + "\n", + "\\begin{equation}\n", + " p(f|\\mathcal{D}) = \\frac{p(\\mathcal{D}|f)p(f)}{p(\\mathcal{D})}\n", + "\\end{equation}\n", + "\n", + "### GP Regression\n", + "\n", + "\\begin{equation}\n", + " y_i = f(\\pmb{x}_i) + \\epsilon_i \\qquad \n", + "\t\\left\\{ \\begin{array}{l}\n", + "\t\tf \\sim \\mathcal{GP}(\\pmb{0}, \\pmb{K}) \\\\\n", + "\t\t\\epsilon_i \\sim \\mathcal{N}(0, \\sigma^2)\n", + "\t\\end{array} \\right.\n", + "\\end{equation}\n", + "\n", + "* Prior $f$ is a GP $\\Leftrightarrow p(\\pmb{f}|\\pmb{X}) = \\mathcal{N}(\\pmb{0}, \\pmb{K})$\n", + "* Likelihood is Gaussian $\\Leftrightarrow p(\\pmb{y}|\\pmb{X},\\pmb{f}) = \\mathcal{N}(\\pmb{f}, \\sigma^2\\pmb{I})$\n", + "* $\\rightarrow p(f|\\mathcal{D})$ is also a GP.\n", + "\n", + "#### Predictive distribution: \n", + "$$p(\\pmb{y}_*|\\pmb{x}_*, \\pmb{X}, \\pmb{y}) = \\int p(\\pmb{y}_{*}| \\pmb{x}_{*}, \\pmb{f}, \\pmb{X}, \\pmb{y}) p(f|\\pmb{X}, \\pmb{y}) d\\pmb{f} = \\mathcal{N}(\\pmb{\\mu}_*, \\pmb{\\Sigma}_*)$$\n", + "\\begin{align}\n", + " \\pmb{\\mu}_* &= \\pmb{K}_{*N} (\\pmb{K}_N + \\sigma^2 \\pmb{I})^{-1} \\pmb{y} \\\\\n", + " \\pmb{\\Sigma}_* &= \\pmb{K}_{**} - \\pmb{K}_{*N} (\\pmb{K}_N + \\sigma^2 \\pmb{I})^{-1} \\pmb{K}_{N*}\n", + "\\end{align}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Learning a GP\n", + "#### Marginal likelihood:\n", + "\\begin{equation}\n", + " p(\\pmb{y}|\\pmb{X}) = \\int p(\\pmb{y}|\\pmb{f},\\pmb{X}) p(\\pmb{f}|\\pmb{X}) d\\pmb{f} = \\mathcal{N}(\\pmb{0}, \\pmb{K} + \\sigma^2\\pmb{I})\n", + "\\end{equation}\n", + "\n", + "By taking the logarithm, and setting $\\pmb{K}_y = (\\pmb{K} + \\sigma^2\\pmb{I})$, we have:\n", + "\n", + "\\begin{equation}\n", + " \\mathcal{L} = \\log p(\\pmb{y}|\\pmb{X}; \\pmb{\\Phi}) = \\underbrace{-\\frac{1}{2} \\pmb{y}^T \\pmb{K}_y^{-1} \\pmb{y}}_{\\mbox{data fit}} \\underbrace{-\\frac{1}{2} \\log |\\pmb{K}_y^{-1}|}_{\\mbox{complexity penalty}} - \\frac{n}{2} \\log 2\\pi\n", + "\\end{equation}\n", + "\n", + "The marginal likelihood (i.e. ML-II) is used to optimize the hyperparameters $\\pmb{\\Phi}$ that defines the covariance function and thus the GP.\n", + "\n", + "\\begin{equation}\n", + " \\pmb{\\Phi}^* = argmax_{\\pmb{\\Phi}} \\log p(\\pmb{y}|\\pmb{X}; \\pmb{\\Phi})\n", + "\\end{equation}\n", + "\n", + "Optimizing the marginal likelihood is more robust than the likelihood as it tries to optimize the complexity of the model, and the fitting of this last one to the observed data." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### GPy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# GP Regression\n", + "# Based on the tutorial: https://github.com/SheffieldML/notebook/blob/master/GPy/GPyCrashCourse.ipynb\n", + "\n", + "# Create dataset\n", + "X = np.random.uniform(-3.0, 3.0, (20,1))\n", + "Y = np.sin(X) + np.random.randn(20,1) * 0.05 \n", + "\n", + "# Create the kernel\n", + "# Reminder 1: The sum of valid kernels gives a valid kernel.\n", + "# Reminder 2: The product of valid kernels gives a valid kernel.\n", + "# Available kernels: RBF, Exponential, Matern32, Matern52, Brownian, Bias, Linear, PeriodicExponential, White.\n", + "kernel = GPy.kern.RBF(input_dim=1, variance=1.0, lengthscale=1.0)\n", + "\n", + "# Create the model\n", + "gp_model = GPy.models.GPRegression(X, Y, kernel)\n", + "\n", + "# Display and plot\n", + "print(\"Before optimization: \", gp_model)\n", + "gp_model.plot()\n", + "plt.show()\n", + "\n", + "# Optimize the model (that is find the 'best' hyperparameters of the kernel matrix)\n", + "# By default, the optimizer is a 2nd order algo: lbfgsb. Others are available such as the scg, ...\n", + "gp_model.optimize(messages=False)\n", + "\n", + "# Display and plot\n", + "print(\"After optimization: \", gp_model)\n", + "gp_model.plot()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Gaussian Process Latent Variable Model (GP-LVM)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# GPLVM\n", + "# Based on the tutorials: \n", + "# http://nbviewer.jupyter.org/github/SheffieldML/notebook/blob/master/GPy/MagnificationFactor.ipynb\n", + "# https://github.com/SheffieldML/notebook/blob/master/lab_classes/gprs/lab4-Copy0.ipynb\n", + "\n", + "# Create dataset\n", + "N = 100\n", + "k1 = GPy.kern.RBF(5, variance=1, lengthscale=1./np.random.dirichlet(np.r_[10,10,10,0.1,0.1]), ARD=True)\n", + "k2 = GPy.kern.RBF(5, variance=1, lengthscale=1./np.random.dirichlet(np.r_[0.1,10,10,10,0.1]), ARD=True)\n", + "X = np.random.normal(0, 1, (N,5))\n", + "A = np.random.multivariate_normal(np.zeros(N), k1.K(X), 10).T\n", + "B = np.random.multivariate_normal(np.zeros(N), k2.K(X), 10).T\n", + "\n", + "Y = np.vstack((A,B))\n", + "\n", + "# latent space dimension\n", + "latent_dim = 2\n", + "\n", + "# Create the kernel\n", + "kernel = GPy.kern.RBF(input_dim=latent_dim, variance=1.0, lengthscale=1.0)\n", + "\n", + "# Create the GPLVM model\n", + "gplvm_model = GPy.models.GPLVM(Y, latent_dim, init='PCA', kernel=kernel)\n", + "\n", + "# Display and plot\n", + "print(\"Before optimization: \", gplvm_model)\n", + "gplvm_model.plot_latent()\n", + "plt.show()\n", + "\n", + "# Optimize the model (that is find the 'best' hyperparameters of the kernel matrix)\n", + "# By default, the optimizer is a 2nd order algo: lbfgsb. Others are available such as the scg, ...\n", + "gplvm_model.optimize(messages=False)\n", + "\n", + "# Display and plot\n", + "print(\"After optimization: \", gplvm_model)\n", + "gplvm_model.plot_latent()\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "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.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/tutorials/machine_learning/README.md b/tutorials/machine_learning/README.md new file mode 100644 index 0000000..dea508f --- /dev/null +++ b/tutorials/machine_learning/README.md @@ -0,0 +1,25 @@ +## Tutorials about machine learning + +This folder contains tutorials about different learning models and algorithms used. + +## References + +Prerequisites: some tutorials require the student to be familiar with multivariate calculus, linear algebra, probability and statistics, information theory, and other mathematical fields. + +Here are references/tutorials that I have watched/read, with their corresponding level: +1. "Machine Learning", by Prof. Andrew Ng. (Easy + practical) +2. "CS188: Introduction to Artificial Intelligence", (Easy + theoretical/practical) +3. "Learning from Data" Yaser (Medium + theoritical) +4. "Machine Learning", by Prof. Nado de Freitas, 2013 (Medium + theoretical) +5. "Deep Learning" (www.deeplearningbook.org), Goodfellow et al., 2016 (Easy) +6. "Gaussian Processes for Machine Learning", Rasmussen and Williams, 2006 (Medium) +7. "Pattern Recognition and Machine Learning", Bishop, 2006 (Hard) +8. "Deep Learning for Computer Vision" (Easy) +9. "Deep Learning for NLP" (Medium) +10. "Reinforcement Learning" (http://www0.cs.ucl.ac.uk/staff/d.silver/web/Teaching.html), Silver, UC London, 2015 (Medium) +11. "An Introduction to Reinforcement Learning", Sutton and Barto, 2018 (Medium) +12. "Deep Reinforcement Learning", (Medium/Hard) + +Are there other stuffs to learn? +Yes, probabilistic graphical models, variational inference, information geometry, etc. + diff --git a/tutorials/math/Linear-Algebra.ipynb b/tutorials/math/Linear-Algebra.ipynb new file mode 100644 index 0000000..d703912 --- /dev/null +++ b/tutorials/math/Linear-Algebra.ipynb @@ -0,0 +1,487 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Linear Algebra" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Vector Space\n", + "\n", + "A ***vector space*** $V$ defined on the field $K$ is ...\n", + "\n", + "It satisfies the 10 following properties:\n", + "1.\n", + "\n", + "A ***vector*** is an element of a vector space." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Linear Transformation\n", + "\n", + "Mapping between 2 vector spaces." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Null and range space\n", + "\n", + "Assume two vector spaces $V$ and $W$, and a linear mapping $\\pmb{A}: V \\rightarrow W$.\n", + "\n", + "* The null (aka kernel) space of a matrix $\\pmb{A}$ is defined as the subspace of $V$ such that:\n", + "\n", + "\\begin{equation}\n", + " \\mathcal{N}(\\pmb{A}) = \\{\\pmb{x} \\in V : \\pmb{Ax} = \\pmb{0}\\} \\subseteq V\n", + "\\end{equation}\n", + "\n", + "* The range (aka column or image) space of a matrix $\\pmb{A}$ is defined as the subspace of $W$ such that:\n", + "\n", + "\\begin{equation}\n", + " \\mathcal{R}(\\pmb{A}) = \\{\\pmb{y} \\in W : \\pmb{Ax} = \\pmb{y}, \\; \\forall \\pmb{x} \\in V \\} \\subseteq W\n", + "\\end{equation}\n", + "\n", + "By the ***Rank–nullity theorem***:\n", + "\\begin{equation}\n", + " dim(\\mathcal{N}(\\pmb{A})) + dim(\\mathcal{R}(\\pmb{A})) = dim(V)\n", + "\\end{equation}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Rank of a matrix" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Transpose" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Determinant\n", + "\n", + "Determinant of a square matrix.\n", + "\n", + "\\begin{equation}\n", + " det : \\mathbb{K}^{N \\times N} \\rightarrow \\mathbb{K}: det(\\pmb{A}) = |\\pmb{A}| \n", + "\\end{equation}\n", + "\n", + "Properties:\n", + "* $det(\\pmb{A}) = det(\\pmb{A}^T)$\n", + "* $det(\\pmb{A}^{-1}) = det(\\pmb{A})^{-1}$\n", + "* If $\\pmb{A}$ and $\\pmb{B}$ are square matrices of same size then $det(\\pmb{A}\\pmb{B}) = det(\\pmb{A}) det(\\pmb{B}) = det(\\pmb{B}) det(\\pmb{A}) = det(\\pmb{B}\\pmb{A})$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Invertible matrix\n", + "\n", + "A square matrix $\\pmb{A} \\in \\mathbb{R}^{N \\times N}$ is invertible, and if $\\pmb{A}\\pmb{A}^{-1} = \\pmb{A}^{-1}\\pmb{A} = \\pmb{I}$.\n", + "\n", + "* $\\pmb{A}$ is invertible.\n", + "* $\\pmb{A}$ is full-rank, i.e. \n", + "* The number 0 is not an eigenvalue of $\\pmb{A}$.\n", + "\n", + "Time complexity: $O(N^3)$\n", + "\n", + "Notes:\n", + "* If you want to know if a matrix $\\pmb{A}$ (which is at least PSD) is PD, you can check if it is invertible. If that's the case, then 0 can not be an eigenvalue of $\\pmb{A}$ and thus the matrix is PD." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Orthogonal Matrices\n", + "\n", + "$\\pmb{A}$ is an orthogonal (square) matrix if $\\pmb{A}^{-1} = \\pmb{A}^T$, and thus $\\pmb{A}^T\\pmb{A} = \\pmb{A}\\pmb{A}^T = \\pmb{I}$. This implies that the columns (resp. rows) of $\\pmb{A}$ are othonormal to each other. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Symmetric Matrices" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Positive (Semi) Definite\n", + "\n", + "* PSD $\\rightarrow$ symmetric.\n", + "\n", + "* All the evals of a PD matrix are positive, and 0 is not one of these.\n", + "* All the evals of a PSD matrix are non-negative." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Norm" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Trace\n", + "\n", + "\\begin{equation}\n", + " tr(\\pmb{A}) = \\sum_i a_{ii}\n", + "\\end{equation}\n", + "\n", + "Properties:\n", + "* The trace is a linear operator: $tr(c_1 \\pmb{A}+ c_2 \\pmb{B}) = c_1 tr(\\pmb{A}) + c_2 tr(\\pmb{B})$\n", + "* The transpose has the same trace: $tr(\\pmb{A}) = tr(\\pmb{A}^T)$\n", + "* Invariance under cyclic permutation: $tr(\\pmb{A}\\pmb{B}\\pmb{C}) = tr(\\pmb{B}\\pmb{C}\\pmb{A}) = tr(\\pmb{C}\\pmb{A}\\pmb{B})$\n", + "* Product: $tr(\\pmb{A}\\pmb{B}) \\neq tr(\\pmb{A})tr(\\pmb{B})$\n", + "* Frobenius Norm: $||\\pmb{A}||_F = \\sqrt{tr(\\pmb{A}^T\\pmb{A})} = \\sqrt{tr(\\pmb{A}\\pmb{A}^T)}$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Matrix Calculus\n", + "\n", + "More information can be found on the [matrix cookbook](https://www.math.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf).\n", + "\n", + "First order:\n", + "\\begin{align}\n", + " \\nabla_{\\pmb{X}} tr(\\pmb{AX}) = \\nabla_{\\pmb{X}} tr(\\pmb{XA}) &= \\pmb{A}^T \\\\\n", + " \\nabla_{\\pmb{X}} tr(\\pmb{AX}^T) = \\nabla_{\\pmb{X}} tr(\\pmb{X}^T\\pmb{A}) &= \\pmb{A}\n", + "\\end{align}\n", + "\n", + "Second order:\n", + "\\begin{align}\n", + " \\nabla_{\\pmb{X}} tr(\\pmb{A}\\pmb{X}\\pmb{B}\\pmb{X}) = \\nabla_{\\pmb{X}} tr(\\pmb{X}\\pmb{B}\\pmb{X}\\pmb{A}) &= (\\pmb{BXA})^T + (\\pmb{AXB})^T \\\\\n", + " \\nabla_{\\pmb{X}} tr(\\pmb{A}\\pmb{X}^T\\pmb{B}\\pmb{X}) = \\nabla_{\\pmb{X}} tr(\\pmb{X}^T\\pmb{B}\\pmb{X}\\pmb{A}) = \\nabla_{\\pmb{X}} tr(\\pmb{B}\\pmb{X}\\pmb{A}\\pmb{X}^T) &= \\pmb{BXA} + (\\pmb{A}\\pmb{X}^T\\pmb{B})^T \\\\\n", + "\\end{align}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Pseudo inverse\n", + "\n", + "This is the generalization of the inverse, in the sense that it can be applied to rectangular matrices. Assume $\\pmb{A} \\in \\mathbb{R}^{N \\times N}$.\n", + "\n", + "Right pseudo-inverse:\n", + "\\begin{equation}\n", + " \\pmb{A}^\\dagger = (\\pmb{A}^T\\pmb{A})^{-1}\\pmb{A}^T\n", + "\\end{equation}\n", + "\n", + "Left pseudo-inverse:\n", + "\\begin{equation}\n", + " ^\\dagger\\pmb{A} = \\pmb{A}^T(\\pmb{A}\\pmb{A}^T)^{-1}\n", + "\\end{equation}\n", + "\n", + "#### Linear Regression (LR)\n", + "\n", + "Assume the input is given by $\\pmb{X} \\in \\mathbb{R}^{N \\times D_x}$, the output by $\\pmb{Y} \\in \\mathbb{R}^{N \\times D_y}$, and the weight matrix by $\\pmb{W} \\in \\mathbb{R}^{D_x \\times D_y}$.\n", + "\n", + "The MSE loss is defined as:\n", + "\n", + "\\begin{equation}\n", + " \\mathcal{L} = ||\\pmb{Y} - \\pmb{XW}||^2_F\n", + "\\end{equation}\n", + "\n", + "Taking the gradient of this loss and setting it to zero gives us the minimum:\n", + "\n", + "\\begin{align}\n", + " \\nabla_{\\pmb{W}} \\mathcal{L} &= \\nabla_{\\pmb{W}} ||\\pmb{Y} - \\pmb{XW}||^2_F \\\\\n", + " &= \\nabla_{\\pmb{W}} tr((\\pmb{Y} - \\pmb{XW})^T(\\pmb{Y} - \\pmb{XW})) \\\\\n", + " &= \\nabla_{\\pmb{W}} [tr(\\pmb{Y}^T\\pmb{Y}) - tr((\\pmb{XW})^T\\pmb{Y}) - tr(\\pmb{Y}^T\\pmb{XW}) + tr(\\pmb{W}^T\\pmb{X}^T\\pmb{X}\\pmb{W})] \\\\\n", + " &= \\nabla_{\\pmb{W}} [tr(\\pmb{Y}^T\\pmb{Y}) - tr(\\pmb{Y}^T\\pmb{XW}) - tr(\\pmb{Y}^T\\pmb{XW}) + tr(\\pmb{W}^T\\pmb{X}^T\\pmb{X}\\pmb{W})] \\\\\n", + " &= -2 \\nabla_{\\pmb{W}} tr(\\pmb{Y}^T\\pmb{XW}) + \\nabla_{\\pmb{W}} tr(\\pmb{W}^T\\pmb{X}^T\\pmb{X}\\pmb{W}) \\\\\n", + " &= -2 (\\pmb{Y}^T\\pmb{X})^T + \\pmb{X}^T\\pmb{X} \\pmb{W} + (\\pmb{X}^T\\pmb{X})^T \\pmb{W} \\\\\n", + " &= -2 \\pmb{X}^T\\pmb{Y} + 2 \\pmb{X}^T\\pmb{X} \\pmb{W} \\\\\n", + " &= 0 \\\\\n", + " \\Leftrightarrow & \\quad (\\pmb{X}^T\\pmb{X}) \\pmb{W} = \\pmb{X}^T\\pmb{Y}\n", + "\\end{align}\n", + "\n", + "If the covariance $(\\pmb{X}^T\\pmb{X})$ is invertible i.e. is PD, then the best set of weights are given by:\n", + "\n", + "\\begin{equation}\n", + " \\pmb{W}^* = (\\pmb{X}^T\\pmb{X})^{-1} \\pmb{X}^T\\pmb{Y} = \\pmb{X}^\\dagger \\pmb{Y} = \\pmb{\\Sigma_{XX}}^{-1}\\pmb{\\Sigma_{XY}}\n", + "\\end{equation}\n", + "\n", + "#### Linear Weighted Regression (LWR)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Covariance\n", + "\n", + "The covariance $\\pmb{C}_{XY}$ is a positive semi-definite (PSD) matrix which captures linear correlation between 2 random variables $X$ and $Y$. The PSD implies that it is symmetric by definition.\n", + "\n", + "* If the cov is invertible, then 0 can not be an eigenvalue of $\\pmb{C}$. This means that it is positive definite (PD)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Eigenvalue and Eigenvectors" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Diagonalizable and Eigendecomposition" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Idempotence\n", + "\n", + "An ***idempotent*** matrix $\\pmb{P}$ is a square 'matrix which, when multiplied by itself, yields itself'.\n", + "\n", + "\\begin{equation}\n", + " \\pmb{P} = \\pmb{PP} = \\pmb{P}^2\n", + "\\end{equation}\n", + "\n", + "Properties:\n", + "* An idempotent matrix (except the identity) is singular (i.e. not full rank).\n", + "* $\\pmb{I} - \\pmb{P}$ is also idempotent.\n", + "* An idempotent matrix is always diagonalizable and its eigenvalues are either 0 or 1.\n", + "* The trace of an idempotent matrix equals the rank of the matrix and thus is always an integer.\n", + "\n", + "In linear regression, the optimal solution of $\\mathcal{L} = ||\\pmb{Y} - \\pmb{XW}||^2_F$ with respect to $\\pmb{W}$ is $\\pmb{W}^* = (\\pmb{X}^T\\pmb{X})^{-1} \\pmb{X}^T\\pmb{Y} = \\pmb{X}^\\dagger \\pmb{Y}$.\n", + "\n", + "The residual error is then given by:\n", + "\n", + "\\begin{equation}\n", + " E = (\\pmb{Y} - \\pmb{XW}^*) = (\\pmb{Y} - \\pmb{X}(\\pmb{X}^T\\pmb{X})^{-1} \\pmb{X}^T\\pmb{Y}) = [\\pmb{I} - \\pmb{X}(\\pmb{X}^T\\pmb{X})^{-1} \\pmb{X}^T] \\pmb{Y} = [\\pmb{I} - \\pmb{X}\\pmb{X}^\\dagger] \\pmb{Y} = \\pmb{Q}\\pmb{Y}\n", + "\\end{equation}\n", + "\n", + "The matrices $\\pmb{P} = \\pmb{X}\\pmb{X}^\\dagger$ and $\\pmb{Q} = [\\pmb{I} - \\pmb{X}\\pmb{X}^\\dagger] = [\\pmb{I} - \\pmb{P}]$ are symmetric and idempotent matrices.\n", + "\n", + "* An idempotent linear operator $\\pmb{P}$ is a projection operator on the range space $\\mathcal{R}(\\pmb{P})$ along its null space $\\mathcal{N}(\\pmb{P})$.\n", + "* $\\pmb{P}$ is an orthogonal projection operator $\\Leftrightarrow$ it is idempotent and symmetric.\n", + "Ex: $\\pmb{P} = \\left[\\begin{array}{cc} \\pmb{I} & \\pmb{0} \\\\ \\pmb{0} & \\pmb{0} \\end{array} \\right]$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Orthogonal Complement and Projection Matrix\n", + "\n", + "\"A ***projection*** is a linear transformation $\\pmb{P}$ from a vector space $V$ to itself such that $\\pmb{P}^2 = \\pmb{P}$ (i.e. $\\pmb{P}$ is idempotent). That is, whenever $\\pmb{P}$ is applied twice to any value, it gives the same result as if it were applied once (idempotent). It leaves its image unchanged.\n", + "\n", + "Let $V$ be a finite dimensional vector space and $\\pmb{P}$ be a projection on $V$. Suppose the subspaces $\\mathcal{R}$ and $\\mathcal{N}$ are the range and null space (aka kernel) of $\\pmb{P}$ respectively. Then $\\pmb{P}$ has the following properties:\n", + "1. $\\pmb{P}$ is idempotent by def (i.e. $\\pmb{P}^2 = \\pmb{P}$)\n", + "2. $\\pmb{P}$ is the identity operator $\\pmb{I}$ on $\\mathcal{R}$ (i.e. $\\forall \\pmb{x} \\in \\mathcal{R}: \\pmb{Px} = \\pmb{x}$)\n", + "3. We have a direct sum $V = \\mathcal{R} \\oplus \\mathcal{N}$. Every vector $\\pmb{x} \\in V$ may be decomposed uniquely as $\\pmb{x} = \\pmb{r} + \\pmb{n}$ with $\\pmb{r} = \\pmb{Px} \\in \\mathcal{R}$ and $\\pmb{n} = \\pmb{x} − \\pmb{Px} = (\\pmb{I} - \\pmb{P}) \\pmb{x} \\in \\mathcal{N}$.\n", + "\n", + "The range and kernel of a projection are complementary, as are $\\pmb{P}$ and $\\pmb{Q} = \\pmb{I} − \\pmb{P}$. The operator $\\pmb{Q}$ is also a projection, and the range and null spaces of $\\pmb{P}$ become the null and range spaces of $\\pmb{Q}$ and vice versa. We say $\\pmb{P}$ is a projection along $\\mathcal{N}$ onto $\\mathcal{R}$, and $\\pmb{Q}$ is a projection along $\\mathcal{R}$ onto $\\mathcal{N}$.\"\n", + "\n", + "##### Orthogonal Projection\n", + "\n", + "When the vector space $V$ has an inner product and is complete (i.e. it is a Hilbert space), the concept of orthogonality can be used. An ***orthogonal projection*** is a projection for which the range $\\mathcal{R}$ and the null space $\\mathcal{N}$ are orthogonal subspaces. That is, $\\forall \\pmb{r} \\in \\mathcal{R}, \\forall \\pmb{n} \\in \\mathcal{N}: \\pmb{r} . \\pmb{n} = 0$.\n", + "\n", + "The ***orthogonal complement*** of a subspace $W$ of a vector space $V$ equipped with a bilinear form $B$ is the set $W^\\perp$ of all vectors in $V$ that are orthogonal to every vector in $W$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### SVD\n", + "\n", + "Any rectangular matrices $\\pmb{A} \\in \\mathbb{R}^{N \\times M}$ can be decomposed into a product of 3 matrices (which can be seen as the building blocks of $\\pmb{A}$):\n", + "\n", + "\\begin{equation}\n", + " \\pmb{A} = \\pmb{U_A}\\pmb{\\Sigma_A}\\pmb{V_A}^T\n", + "\\end{equation}\n", + "where:\n", + "* $\\pmb{U_A} \\in \\mathbb{R}^{N \\times N}$ is an orthogonal matrix (i.e. $\\pmb{U_A}^T = \\pmb{U_A}^{-1}$ thus $\\pmb{U_A}\\pmb{U_A}^T = \\pmb{U_A}^T\\pmb{U_A} = \\pmb{I}$). The columns of $\\pmb{U_A}$ contains the eigenvectors of the PSD (thus symmetric) $\\pmb{A}\\pmb{A}^T$, and are also known as the left-singular vectors of $\\pmb{A}$. Its first columns associated with non-zero singular values span the range of $\\pmb{A}$. Note that because $\\pmb{U_A}$ is an orthogonal matrix, it is invertible and thus has full rank, i.e. its columns form a basis and span $\\mathbb{R}^N$.\n", + "* $\\pmb{\\Sigma_A} \\in \\mathbb{R}^{N \\times M}$ is a rectangular matrix. The upper left submatrix is a diagonal matrix of size $r \\times r$ with $r = \\min(N,M)$ while the rest is filled with zeros. The diagonal elements are the singular values ordered by descending order. The singular values (SVs) $\\sigma_i$ are equals to the square roots of the eigenvalues $\\lambda_i$ of $\\pmb{A}\\pmb{A}^T$ and $\\pmb{A}^T\\pmb{A}$. The rank of $\\pmb{A}$ is given by the number of SVs different from 0.\n", + "* $\\pmb{V_A} \\in \\mathbb{R}^{M \\times M}$ is an orthogonal matrix (i.e. $\\pmb{V_A}^T = \\pmb{V_A}^{-1}$ thus $\\pmb{V_A}\\pmb{V_A}^T = \\pmb{V_A}^T\\pmb{V_A} = \\pmb{I}$). The columns of $\\pmb{V_A}$ contains the eigenvectors of the PSD (thus symmetric) $\\pmb{A}^T\\pmb{A}$, and are also known as the right-singular vectors of $\\pmb{A}$. Its last columns associated with vanishing singular values $\\sigma_i = 0$ span the null space of $\\pmb{A}$. Note that because $\\pmb{V_A}$ is an orthogonal matrix, it is invertible and thus has full rank, i.e. its columns form a basis and span $\\mathbb{R}^M$.\n", + "\n", + "Few notes:\n", + "* SVD can be seen as a generalization of eigendecomposition.\n", + "* SVD can be compressed such that $\\pmb{A} = \\pmb{\\tilde{U}_A}\\pmb{\\tilde{\\Sigma}_A}\\pmb{\\tilde{V_A}^T}$ with $\\pmb{\\tilde{U}_A} \\in \\mathbb{R}^{N \\times r}$, $\\pmb{\\tilde{\\Sigma}_A} \\in \\mathbb{R}^{r \\times r}$, $\\pmb{V_A} \\in \\mathbb{R}^{M \\times r}$, and $r=\\min(N,M)$.\n", + "* \n", + "\n", + "Some properties:\n", + "* **Existence**:\n", + "* **Uniqueness**:\n", + "* **Transpose**: $\\pmb{A}^T = (\\pmb{U_A}\\pmb{\\Sigma_A}\\pmb{V_A}^T)^T = \\pmb{V_A}\\pmb{\\Sigma_A}^T\\pmb{U_A}^T $\n", + "* If the matrix A is symmetric, then it has real eigenvalues.\n", + "* **Pseudo-inverse**: The pseudo-inverse of $\\pmb{A}$ is given by $\\pmb{A^\\dagger} = \\pmb{V_A}\\pmb{\\Sigma_A}^{-1}\\pmb{U_A}^T$ with $\\pmb{\\Sigma_A}^{-1}$ containing the inverse of singular values on its diagonal.\n", + "* **Null and range space** of $\\pmb{A}$: $\\mathcal{N}(\\pmb{A}) \\equiv \\mathcal{R}^\\perp(\\pmb{A}^T)$ and $\\mathcal{R}(\\pmb{A}) \\equiv \\mathcal{N}^\\perp(\\pmb{A}^T)$. This can be seen by applying SVD on $\\pmb{A}$ and $\\pmb{A}^T$. The last columns of $\\pmb{V_A}$ associated with vanishing SVs spans the null-space of $\\pmb{A}$, and the first columns of $\\pmb{V_A}$ associated with non-zero SVs spans the range space of $\\pmb{A}^T$. Because of the orthogonality property of $\\pmb{V_A}$, we have the $r$ first columns are orthogonals to the $M-r$ last columns. $dim(\\mathcal{R}(\\pmb{A})) + dim(\\mathcal{N}(\\pmb{A})) = M$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Block Matrices" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### PCA\n", + "\n", + "The (MSE) loss minimized by PCA is:\n", + "\n", + "\\begin{equation}\n", + " \\mathcal{L} = ||\\pmb{X} - \\pmb{XWW}^T||^2_F\n", + "\\end{equation}\n", + "\n", + "where is the $\\pmb{\\tilde{X}} = \\pmb{XW}$ is the projected data on the lower dimensional space, and $\\pmb{\\hat{X}} = \\pmb{\\tilde{X}W}^T$ is the data projected back to the original space.\n", + "\n", + "Taking the gradient of this loss with respect to $\\pmb{W}$ and setting it to $0$ gives us the minimum:\n", + "\n", + "\\begin{align}\n", + " \\nabla_{\\pmb{W}} \\mathcal{L} &= \\nabla_{\\pmb{W}} tr((\\pmb{X} - \\pmb{XWW}^T)^T(\\pmb{X} - \\pmb{XWW}^T)) \\\\\n", + " &= \\nabla_{\\pmb{W}} tr((\\pmb{X}^T - \\pmb{WW}^T \\pmb{X}^T) (\\pmb{X} - \\pmb{XWW}^T)) \\\\\n", + " &= \\nabla_{\\pmb{W}} [tr(\\pmb{X}^T\\pmb{X}) - tr(\\pmb{WW}^T \\pmb{X}^T\\pmb{X}) - tr(\\pmb{X}^T\\pmb{X}\\pmb{WW}^T) + tr(\\pmb{WW}^T\\pmb{X}^T\\pmb{XWW}^T)] \\\\\n", + " &= \\nabla_{\\pmb{W}} [- 2 tr(\\pmb{W}^T \\pmb{\\Sigma_{XX}} \\pmb{W}) + tr(\\pmb{WW}^T \\pmb{\\Sigma_{XX}} \\pmb{WW}^T)] \\\\\n", + " &= - 2 (\\pmb{\\Sigma_{XX}} \\pmb{W} + \\pmb{\\Sigma_{XX}}^T \\pmb{W}) + (\\pmb{W}^T \\pmb{\\Sigma_{XX}} \\pmb{WW}^T)^T + (\\pmb{\\Sigma_{XX}} \\pmb{WW}^T)\\pmb{W} + (\\pmb{WW}^T \\pmb{\\Sigma_{XX}})^T (\\pmb{W}^T)^T + (\\pmb{WW}^T \\pmb{\\Sigma_{XX}} \\pmb{W}) \\\\\n", + " &= - 4 \\pmb{\\Sigma_{XX}} \\pmb{W} + \\pmb{WW}^T \\pmb{\\Sigma_{XX}}^T \\pmb{W} + \\pmb{\\Sigma_{XX}} \\pmb{WW}^T\\pmb{W} + \\pmb{\\Sigma_{XX}}^T\\pmb{WW}^T \\pmb{W} + \\pmb{WW}^T \\pmb{\\Sigma_{XX}} \\pmb{W} \\\\\n", + " &= - 4 \\pmb{\\Sigma_{XX}} \\pmb{W} + 2 \\pmb{WW}^T \\pmb{\\Sigma_{XX}} \\pmb{W} + 2 \\pmb{\\Sigma_{XX}} \\pmb{WW}^T\\pmb{W} \\\\\n", + " &= 0 \\\\\n", + " \\Leftrightarrow & \\quad \\pmb{WW}^T \\pmb{\\Sigma_{XX}} \\pmb{W} + \\pmb{\\Sigma_{XX}} \\pmb{WW}^T\\pmb{W} = 2 \\pmb{\\Sigma_{XX}} \\pmb{W} \\\\\n", + " \\Leftrightarrow & \\quad \\pmb{WW}^T \\pmb{\\Sigma_{XX}} + \\pmb{\\Sigma_{XX}} \\pmb{WW}^T = 2 \\pmb{\\Sigma_{XX}} \\\\\n", + " \\Leftrightarrow & \\quad \\pmb{WW}^T + \\pmb{\\Sigma_{XX}} \\pmb{WW}^T \\pmb{\\Sigma_{XX}}^{-1} = 2 \\pmb{I} \\\\\n", + " \\Leftrightarrow & \\quad \\pmb{\\Sigma_{XX}} = (2 \\pmb{I} - \\pmb{WW}^T) \\pmb{\\Sigma_{XX}} (\\pmb{WW}^T)^{-1} \\\\\n", + " \\Leftrightarrow & \\quad \\pmb{Q\\Lambda Q}^T = (2 \\pmb{I} - \\pmb{WW}^T) \\pmb{Q\\Lambda Q}^T (\\pmb{WW}^T)^{-1} \\\\\n", + " \\Leftrightarrow & \\quad \\left\\{ \\begin{array}{l} (2 \\pmb{I} - \\pmb{WW}^T) \\pmb{Q} = \\pmb{Q} \\\\ \\pmb{Q}^T (\\pmb{WW}^T)^{-1} = \\pmb{Q}^T \\: \\Leftrightarrow \\: \\pmb{QQ}^T = \\pmb{WW}^T \\end{array} \\right. \\\\\n", + "\\end{align}\n", + "\n", + "Algo using the covariance:\n", + "1. subtract the mean $\\pmb{X}$\n", + "2. compute the covariance matrix $\\pmb{\\Sigma_{XX}} = \\pmb{X}^T\\pmb{X}$\n", + "3. compute the eigendecomposition of $\\pmb{\\Sigma_{XX}}$, i.e. compute $\\pmb{Q}$ and $\\pmb{\\Lambda}$ such that $\\pmb{\\Sigma_{XX}} = \\pmb{Q} \\pmb{\\Lambda} \\pmb{Q}^T$.\n", + "4. return the sorted evals and the corresponding evecs\n", + "\n", + "Algo using SVD:\n", + "1. substract the mean from $\\pmb{X}$\n", + "2. compute SVD of $\\pmb{X}$, i.e. $\\pmb{X} = \\pmb{U_X \\Sigma_X V_X}^T$. The eigenvectors are given by $\\pmb{Q} = \\pmb{V_X}$ and the evals by $\\pmb{\\Lambda} = \\pmb{\\Sigma_X}^2$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Recursive PCA/SVD" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Hierarchical PCA/SVD" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Linear FNN with 1 hidden layer\n", + "\n", + "Assume a linear Feedforward Neural Network (FNN) with 1 input, 1 hidden, and 1 output layer. Assume the input data is given by $\\pmb{X} \\in \\mathbb{R}^{N \\times D_x}$, the output data by $\\pmb{Y} \\in \\mathbb{R}^{N \\times D_y}$, the hidden data by $\\pmb{H} \\in \\mathbb{R}^{N \\times D_h}$, the input-hidden weight matrix by $\\pmb{W_1} \\in \\mathbb{R}^{D_x \\times D_h}$ and the hidden-output weight matrix $\\pmb{W_2} \\in \\mathbb{R}^{D_h \\times D_y}$. These variables are related by the following relationships:\n", + "\n", + "\\begin{equation}\n", + " \\pmb{H} = \\pmb{X} \\pmb{W_1} \\qquad \\mbox{and} \\qquad \\pmb{Y} = \\pmb{H} \\pmb{W_2}\n", + "\\end{equation}\n", + "\n", + "The MSE loss is thus defined as:\n", + "\n", + "\\begin{equation}\n", + " \\mathcal{L} = ||\\pmb{Y} - \\pmb{XW_1W_2}||^2_F\n", + "\\end{equation}\n", + "\n", + "Taking the gradient of this loss with respect to $\\pmb{W_1}$ and $\\pmb{W_2}$, and setting these to zero gives us the minimum:\n", + "\n", + "\\begin{align}\n", + " \\nabla_{\\pmb{W_1}} \\mathcal{L} &= \\nabla_{\\pmb{W_1}} ||\\pmb{Y} - \\pmb{XW_1W_2}||^2_F \\\\\n", + " &= \\nabla_{\\pmb{W_1}} tr((\\pmb{Y} - \\pmb{XW_1W_2})^T(\\pmb{Y} - \\pmb{XW_1W_2})) \\\\\n", + " &= \\nabla_{\\pmb{W_1}} [tr(\\pmb{Y}^T\\pmb{Y}) - tr((\\pmb{XW_1W_2})^T\\pmb{Y}) - tr(\\pmb{Y}^T\\pmb{XW_1W_2}) + tr(\\pmb{W}^T\\pmb{X}^T\\pmb{X}\\pmb{W})] \\\\\n", + " &= \\nabla_{\\pmb{W_1}} [tr(\\pmb{Y}^T\\pmb{Y}) - tr(\\pmb{Y}^T\\pmb{XW_1W_2}) - tr(\\pmb{Y}^T\\pmb{XW_1W_2}) + tr(\\pmb{(W_1W_2)}^T\\pmb{X}^T\\pmb{X}\\pmb{(W_1W_2)})] \\\\\n", + " &= \\nabla_{\\pmb{W_1}} [-2 tr(\\pmb{Y}^T\\pmb{XW_1W_2}) + tr(\\pmb{W_2}^T\\pmb{W_1}^T\\pmb{X}^T\\pmb{X}\\pmb{W_1}\\pmb{W_2})] \\\\\n", + " &= \\nabla_{\\pmb{W_1}} [-2 tr(\\pmb{\\Sigma_{YX}}\\pmb{W_1W_2}) + tr(\\pmb{W_2}^T\\pmb{W_1}^T\\pmb{\\Sigma_{XX}}\\pmb{W_1}\\pmb{W_2})] \\\\\n", + " &= -2 (\\pmb{W_2}\\pmb{\\Sigma_{YX}})^T + \\pmb{\\Sigma_{XX}}\\pmb{W_1W_2W_2}^T + \\pmb{\\Sigma_{XX}}^T\\pmb{W_1}(\\pmb{W_2W_2}^T)^T \\\\\n", + " &= -2 \\pmb{\\Sigma_{XY}} \\pmb{W_2}^T + 2 \\pmb{\\Sigma_{XX}}\\pmb{W_1W_2W_2}^T \\\\\n", + " &= 0 \\\\\n", + " \\Leftrightarrow & \\quad \\pmb{\\Sigma_{XX}}\\pmb{W_1W_2W_2}^T = \\pmb{\\Sigma_{XY}} \\pmb{W_2}^T \\\\\n", + " \\Leftrightarrow & \\quad \\pmb{W_1} = \\pmb{\\Sigma_{XX}}^{-1} \\pmb{\\Sigma_{XY}} \\pmb{W_2}^T (\\pmb{W_2W_2}^T)^{-1}\n", + "\\end{align}\n", + "\n", + "and \n", + "\n", + "\\begin{align}\n", + " \\nabla_{\\pmb{W_2}} \\mathcal{L} &= \\nabla_{\\pmb{W_2}} [-2 tr(\\pmb{\\Sigma_{YX}}\\pmb{W_1W_2}) + tr(\\pmb{W_2}^T\\pmb{W_1}^T\\pmb{\\Sigma_{XX}}\\pmb{W_1}\\pmb{W_2})] \\\\\n", + " &= -2 (\\pmb{\\Sigma_{YX}}\\pmb{W_1})^T + \\pmb{W_1}^T\\pmb{\\Sigma_{XX}}\\pmb{W_1}\\pmb{W_2} + (\\pmb{W_1}^T\\pmb{\\Sigma_{XX}}\\pmb{W_1})^T\\pmb{W_2} \\\\\n", + " &= -2 \\pmb{W_1}^T\\pmb{\\Sigma_{XY}} + 2 \\pmb{W_1}^T\\pmb{\\Sigma_{XX}}\\pmb{W_1}\\pmb{W_2} \\\\\n", + " &= 0 \\\\\n", + " \\Leftrightarrow & \\quad \\pmb{W_1}^T\\pmb{\\Sigma_{XX}}\\pmb{W_1}\\pmb{W_2} = \\pmb{W_1}^T \\pmb{\\Sigma_{XY}} \\\\\n", + " \\Leftrightarrow & \\quad \\pmb{W_2} = (\\pmb{W_1}^T\\pmb{\\Sigma_{XX}}\\pmb{W_1})^{-1} \\pmb{W_1}^T \\pmb{\\Sigma_{XY}} \\\\\n", + "\\end{align}\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "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.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/tutorials/math/PCA.ipynb b/tutorials/math/PCA.ipynb new file mode 100644 index 0000000..9a2f230 --- /dev/null +++ b/tutorials/math/PCA.ipynb @@ -0,0 +1,1144 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# PCA\n", + "\n", + "In this notebook, I will review a little bit about [PCA](https://arxiv.org/abs/1404.1100), implement [recursive PCA](http://www.sciencedirect.com/science/article/pii/S0959152400000226), how PCA can be viewed as an optimization problem, and implement a constrained version of this optimization process for PCA applied on time-series by including a roughness penalty.\n", + "\n", + "[Generalized PCA](https://arxiv.org/abs/1202.4002) will not be considered after careful consideration. Note: \"Generalized PCA\" is \"an algebro-geometric solution to the problem of segmenting an unknown number of subspaces of unknown and varying dimensions from sample data points\"." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from mpl_toolkits.mplot3d import Axes3D\n", + "from matplotlib.patches import FancyArrowPatch\n", + "from mpl_toolkits.mplot3d import proj3d\n", + "\n", + "from sklearn.decomposition import PCA\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class Arrow3D(FancyArrowPatch):\n", + " def __init__(self, xs, ys, zs, *args, **kwargs):\n", + " FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)\n", + " self._verts3d = xs, ys, zs\n", + "\n", + " def draw(self, renderer):\n", + " xs3d, ys3d, zs3d = self._verts3d\n", + " xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)\n", + " self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))\n", + " FancyArrowPatch.draw(self, renderer)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's first start to apply PCA on the following toy example." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Toy data\n", + "X = np.array([[-1, -1],\n", + " [-2, -1],\n", + " [-3, -2],\n", + " [1, 1],\n", + " [2, 1],\n", + " [3, 2]], dtype=np.float64) # NxT\n", + "\n", + "# Plot\n", + "plt.figure(figsize=(5,4))\n", + "plt.scatter(X[:,0], X[:,1], color='b')\n", + "plt.arrow(0, 0, 1, 0, length_includes_head = True, head_width = 0.15, color='k')\n", + "plt.arrow(0, 0, 0, 1, length_includes_head = True, head_width = 0.15, color='')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "### PCA\n", + "n_components = 2 # number of principal axes that we want to keep\n", + "\n", + "## Let's compute PCA from scratch\n", + "# 1. Center the data\n", + "mean = X.mean(axis=0)\n", + "X -= mean # note that the data was already centered\n", + "N = X.shape[0]\n", + "\n", + "# 2. Compute the covariance matrix\n", + "CovX = 1./(N-1) * X.T.dot(X) # TxT (same as np.cov(X, rowvar=False)))\n", + "\n", + "# 3. Compute the eigenvectors of this covariance matrix\n", + "# np.linalg.eigh is more efficient than np.linalg.eig for symmetric matrix\n", + "evals, evecs = np.linalg.eigh(CovX)\n", + "\n", + "# 4. Sort the eigenvalues (in decreasing order) and eigenvectors\n", + "idx = np.argsort(evals)[::-1]\n", + "evals = evals[idx]\n", + "evecs = evecs[:,idx]\n", + "\n", + "# 5. Form the projection matrix\n", + "P = evecs[:,:n_components]\n", + "print(P)\n", + "\n", + "# 5. Project the data using the projection matrix\n", + "# This is the same as rotating the matrix X using P\n", + "Y = X.dot(P)\n", + " \n", + "# 6. Compare it with standard PCA\n", + "pca = PCA(n_components=n_components)\n", + "pca = pca.fit(X)\n", + "Xnew = pca.transform(X)\n", + "\n", + "print(pca.components_.T)\n", + "print(np.allclose(Xnew, Y))\n", + "\n", + "# 7. Plot the data\n", + "plt.figure(figsize=(10,4))\n", + "plt.subplot(1,2,1)\n", + "plt.title('PCA: eigenvalues')\n", + "plt.bar(np.array([0.,0.1]), evals, width=0.1)\n", + "plt.xlim(0.,1.)\n", + "\n", + "plt.subplot(1,2,2)\n", + "plt.title('PCA: data and eigenvectors')\n", + "plt.scatter(X[:,0], X[:,1], color='b')\n", + "plt.arrow(0, 0, np.sqrt(evals[0])*P[0,0], np.sqrt(evals[0])*P[1,0],\n", + " length_includes_head = True, head_width = 0.15, color='b')\n", + "plt.arrow(0, 0, np.sqrt(evals[1])*P[0,1], np.sqrt(evals[1])*P[1,1],\n", + " length_includes_head = True, head_width = 0.15, color='b')\n", + "plt.scatter(Y[:,0], Y[:,1], color='r')\n", + "plt.arrow(0, 0, np.sqrt(evals[0]), 0, length_includes_head = True, head_width = 0.15, color='r')\n", + "plt.arrow(0, 0, 0, np.sqrt(evals[1]), length_includes_head = True, head_width = 0.15, color='r')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The covariance matrix can be recovered from these eigenvalues and eigenvectors." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Using the eigenvectors and eigenvalues, we can of course recover the covariance matrix\n", + "# Note: if P is not square, we have to fill it up.\n", + "np.allclose(CovX, P.dot(np.diag(evals)).dot(P.T))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Few properties:\n", + "\n", + "* Applying PCA several times is the same as applying it one time. This is because PCA diagonalizes our matrix, thus applying PCA on a diagonal matrix will result in the same matrix.\n", + "* Applying PCA on a part of the data and another PCA on the other part, then applying PCA on the concatenation of both do not result in the same matrix as applying PCA on the whole data.\n", + "* Applying PCA on a rotated matrix does not give the same result as applying PCA on the initial matrix." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Here is the general method\n", + "\n", + "def pca(X, normalize=False, copy=True):\n", + " if copy:\n", + " X = np.copy(X)\n", + "\n", + " # 1. Center the data\n", + " mean = X.mean(axis=0)\n", + " X -= mean\n", + " N = X.shape[0]\n", + " \n", + " if normalize:\n", + " X /= X.std(axis=0)\n", + "\n", + " # 2. Compute the covariance matrix\n", + " CovX = 1./(N-1) * X.T.dot(X) # TxT (same as np.cov(X, rowvar=False)))\n", + "\n", + " # 3. Compute the eigenvectors of this covariance matrix\n", + " # np.linalg.eigh is more efficient than np.linalg.eig for symmetric matrix\n", + " evals, evecs = np.linalg.eigh(CovX)\n", + "\n", + " # 4. Sort the eigenvalues (in decreasing order) and eigenvectors\n", + " idx = np.argsort(evals)[::-1]\n", + " evals, evecs = evals[idx], evecs[:,idx]\n", + "\n", + " return evals, evecs" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Applying PCA on time series\n", + "\n", + "A fundamental question when applying PCA on time series is how to visualize this high dimensional data. Indeed, a sample $\\pmb{x}(t) \\in \\mathbb{R}^T$. One way is to plot this $\\pmb{x}(t)$ where the x-axis is the time, and y-axis is $x(t)$. Each time $t_i$ $(\\forall i \\in \\{0,...,T\\})$ represents a dimension. By plotting a vertical line at time $t=t_i$, we can see the variance in this particular dimension.\n", + "\n", + "For an infinite vector, or function, check about *Functional PCA*." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Recursive PCA\n", + "\n", + "Let's now apply **recursive PCA** on this toy example, with 3 new samples coming at different time steps." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Let's augment our matrix with 3 new samples\n", + "Xs = np.array([[-1,1],\n", + " [-3,0],\n", + " [-4,-1]], dtype=np.float64)\n", + "\n", + "n_components = 2\n", + "k = float(N)\n", + "R = CovX\n", + "X_aug = X\n", + "mean = X.mean(axis=0).reshape(-1,1)\n", + "print(evals)\n", + "for x in Xs:\n", + " x = x.reshape(-1,1)\n", + " X_aug = np.vstack((X_aug, x.T))\n", + " X_aug1 = X_aug - X_aug.mean(axis=0)\n", + " pca = PCA(n_components=n_components)\n", + " pca = pca.fit(X_aug1)\n", + " #print(pca.get_covariance())\n", + "\n", + " # Recursive PCA\n", + " new_mean = k/(k+1) * mean + 1./(k+1) * x\n", + " diff_mean = (new_mean - mean)\n", + " R = (k-1)/k * R + diff_mean.dot(diff_mean.T) + 1./k * (x-new_mean).dot((x-new_mean).T)\n", + " #print(R)\n", + " print(\"The cov of the whole augmented matrix is equal to the recursive cov: {0}\".format(\n", + " np.allclose(pca.get_covariance(), R)))\n", + " k+=1\n", + " mean = new_mean\n", + " \n", + " evals = np.linalg.eigh(R)[0]\n", + " idx = np.argsort(evals)[::-1]\n", + " evals = evals[idx]\n", + " print(evals)\n", + "\n", + "# Compute the new projection matrix\n", + "evals, evecs = np.linalg.eigh(R)\n", + "idx = np.argsort(evals)[::-1]\n", + "evals, evecs = evals[idx], evecs[:,idx]\n", + "P = evecs[:,:n_components]\n", + "Y = X_aug.dot(P)\n", + " \n", + "# Plot the data\n", + "plt.figure(figsize=(10,4))\n", + "plt.subplot(1,2,1)\n", + "plt.title('PCA: eigenvalues')\n", + "plt.bar(np.array([0.,0.1]), evals, width=0.1)\n", + "plt.xlim(0.,1.)\n", + "\n", + "plt.subplot(1,2,2)\n", + "plt.title('PCA: data and eigenvectors')\n", + "plt.scatter(X_aug[:,0], X_aug[:,1], color='b')\n", + "plt.arrow(0, 0, np.sqrt(evals[0])*P[0,0], np.sqrt(evals[0])*P[1,0],\n", + " length_includes_head = True, head_width = 0.15, color='b')\n", + "plt.arrow(0, 0, np.sqrt(evals[1])*P[0,1], np.sqrt(evals[1])*P[1,1],\n", + " length_includes_head = True, head_width = 0.15, color='b')\n", + "plt.scatter(Y[:,0], Y[:,1], color='r')\n", + "plt.arrow(0, 0, np.sqrt(evals[0]), 0, length_includes_head = True, head_width = 0.15, color='r')\n", + "plt.arrow(0, 0, 0, np.sqrt(evals[1]), length_includes_head = True, head_width = 0.15, color='r')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's now add a sample that can not be modeled by a linear combination of the principal axes, i.e. which is orthogonal to the current covariance matrix. Then, as usual, let's apply recursive PCA on it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Let's add a sample that can not be modeled by a linear combination of the principal axes\n", + "# i.e. which is orthogonal to the current covariance matrix.\n", + "# Then, let's apply recursive PCA on it.\n", + "\n", + "# New 3D sample\n", + "x = np.array([1,-1,1], dtype=np.float64).reshape(-1,1)\n", + "\n", + "# Reshaping previous values (pad a column/row of zeros)\n", + "X_aug = np.pad(X_aug, ((0,0), (0,1)), mode='constant') # Nx(T+1)\n", + "mean = np.pad(mean, ((0,1),(0,0)), mode='constant')\n", + "R = np.pad(R, ((0,1),(0,1)), mode='constant')\n", + "\n", + "# Adding new sample and compute mean\n", + "X_aug = np.vstack((X_aug, x.T))\n", + "X_aug1 = X_aug - X_aug.mean(axis=0)\n", + "\n", + "# Use sklearn PCA\n", + "n_components = 3\n", + "pca = PCA(n_components=n_components)\n", + "pca = pca.fit(X_aug1)\n", + "#print(pca.get_covariance())\n", + "\n", + "# Recursive PCA\n", + "new_mean = k/(k+1) * mean + 1./(k+1) * x\n", + "diff_mean = (new_mean - mean)\n", + "R = (k-1)/k * R + diff_mean.dot(diff_mean.T) + 1./k * (x-new_mean).dot((x-new_mean).T)\n", + "#print(R)\n", + "print(\"The cov of the whole augmented matrix is equal to the recursive cov: {0}\".format(\n", + " np.allclose(pca.get_covariance(), R)))\n", + "k+=1\n", + "mean = new_mean\n", + "#print('-'*30)\n", + "\n", + "# Compute the new projection matrix\n", + "evals, evecs = np.linalg.eigh(R)\n", + "idx = np.argsort(evals)[::-1]\n", + "evals, evecs = evals[idx], evecs[:,idx]\n", + "P = evecs[:,:n_components]\n", + "Y = X_aug.dot(P)\n", + " \n", + "# Plot the data\n", + "fig = plt.figure(figsize=(10,4))\n", + "plt.subplot(1,2,1)\n", + "plt.title('PCA: eigenvalues')\n", + "plt.bar(np.array([0.,0.1,0.2]), evals, width=0.1)\n", + "plt.xlim(0.,1.)\n", + "\n", + "ax = fig.add_subplot(122, projection='3d')\n", + "ax.set_title('PCA: data and eigenvectors')\n", + "ax.scatter(X_aug[:,0], X_aug[:,1], X_aug[:,2])\n", + "# From https://stackoverflow.com/questions/22867620/putting-arrowheads-on-vectors-in-matplotlibs-3d-plot\n", + "for v in evecs:\n", + " a = Arrow3D([0., v[0]], [0., v[1]], [0., v[2]],\n", + " mutation_scale=20, lw=1, arrowstyle=\"-|>\", color=\"b\")\n", + " ax.add_artist(a)\n", + "ax.scatter(Y[:,0], Y[:,1], Y[:,2], color='r')\n", + "a = Arrow3D([0., evals[0]], [0., evals[1]], [0., evals[2]],\n", + " mutation_scale=20, lw=1, arrowstyle=\"-|>\", color=\"r\")\n", + "ax.add_artist(a)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "## PCA as an Optimization Problem\n", + "\n", + "PCA can be viewed as an optimization problem in 2 different ways. Theses 2 approaches are equivalent.\n", + "1. Maximize the variance of the projected data.\n", + "2. Minimize the reconstruction error in a least-square sense.\n", + "\n", + "Mathematically, here is the optimization problem that we are trying to solve:\n", + "\n", + "\\begin{equation}\n", + " \\max_{\\pmb{v_i}} \\: \\pmb{v_i}^T \\pmb{X}^T \\pmb{X v_i} \\quad \\mbox{subj. to} \\quad \\begin{array}{l} \\pmb{v_i}^T \\pmb{v_i} = 1 \\\\ \\pmb{v_i}^T \\pmb{v_j} = 0\n", + "\\end{array},\n", + "\\end{equation}\n", + "$\\forall i \\in \\{1,...,D\\}, \\forall 1\\leq j < i$.\n", + "\n", + "Nice references:\n", + "* [What is the objective fct of PCA? (StackExchange)](https://stats.stackexchange.com/questions/10251/what-is-the-objective-function-of-pca)\n", + "* [PCA objective function: what is the connection between maximizing variance and minimizing error? (StackExchange)](https://stats.stackexchange.com/questions/32174/pca-objective-function-what-is-the-connection-between-maximizing-variance-and-m)\n", + "* [Everything you did and didn't know about PCA (blog)](http://alexhwilliams.info/itsneuronalblog/2016/03/27/pca/)\n", + "* [\"PCA and Optimization - A Tutorial\", 2015 (paper)](http://scholarscompass.vcu.edu/cgi/viewcontent.cgi?article=1006&context=ssor_pubs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# data\n", + "# Toy data\n", + "X = np.array([[-1, -1],\n", + " [-2, -1],\n", + " [-3, -2],\n", + " [1, 1],\n", + " [2, 1],\n", + " [3, 2]], dtype=np.float64) # NxT\n", + "N = X.shape[0]\n", + "\n", + "# PCA\n", + "evals, evecs = pca(X)\n", + "print(evals)\n", + "print(evecs)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Using Scipy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.optimize import minimize\n", + "\n", + "# PCA as an optimization\n", + "\n", + "# cache the 'covariance' matrix\n", + "C = X.T.dot(X)/(N-1)\n", + "\n", + "# define objective function to MINIMIZE\n", + "f = lambda u: -(u.T.dot(C)).dot(u)\n", + "\n", + "# define initial guess\n", + "x0 = np.zeros((2,1))\n", + "\n", + "# define optimization method\n", + "# By default, it will be 'BFGS', 'L-BFGS-B', or 'SLSQP' depending on the constraints and bounds\n", + "method = None\n", + "\n", + "# define constraints\n", + "constraints = [{'type': 'eq', 'fun': lambda u: u.T.dot(u) - 1}]\n", + "\n", + "# Minimize --> it returns an instance of OptimizeResult\n", + "u1 = minimize(f, x0, method=method, constraints=constraints)\n", + "#print(u1)\n", + "u1 = u1.x.reshape(-1,1) # get solution\n", + "\n", + "# Add constraint\n", + "constraints.append({'type': 'eq', 'fun': lambda u: u1.T.dot(u)})\n", + "u2 = minimize(f, x0, method=method, constraints=constraints)\n", + "#print(u2)\n", + "u2 = u2.x.reshape(-1,1)\n", + "\n", + "# stack the optimized vector found\n", + "P = np.hstack((u1,u2))\n", + "print(P)\n", + "\n", + "# Plot\n", + "plt.title('PCA: data and eigenvectors')\n", + "plt.scatter(X[:,0], X[:,1], color='b')\n", + "plt.arrow(0, 0, P[0,0], P[1,0], length_includes_head = True, head_width = 0.15, color='b')\n", + "plt.arrow(0, 0, P[0,1], P[1,1], length_includes_head = True, head_width = 0.15, color='g')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define PCA optimization method\n", + "def pca_scipy(X, rough_param=0.0, normalize=False, copy=True):\n", + " \"\"\"\n", + " Compute PCA on the given data using an optimization process.\n", + " \"\"\"\n", + " if copy:\n", + " X = np.copy(X)\n", + "\n", + " # center the data\n", + " mean = X.mean(axis=0)\n", + " X -= mean\n", + " N = X.shape[0]\n", + " T = X.shape[1]\n", + "\n", + " # normalize\n", + " if normalize:\n", + " X /= X.std(axis=0)\n", + "\n", + " \n", + " # compute 'covariance' matrix and cache it\n", + " N = X.shape[0]\n", + " C = X.T.dot(X)/(N-1)\n", + "\n", + " # define objective function to MINIMIZE\n", + " #f = lambda u: -(u.T.dot(C)).dot(u)\n", + " def f(u):\n", + " pen = 0\n", + " if rough_param != 0 and u.size > 2:\n", + " ddu = np.diff(np.diff(u))\n", + " rough_pen = (ddu**2).sum()\n", + " pen = rough_param*rough_pen\n", + " #if u.size > 4:\n", + " # ddddu = np.diff(np.diff(ddu))\n", + " # rough_pen = (ddddu**2).sum()\n", + " # pen += rough_param*rough_pen\n", + " loss = -(u.T.dot(C)).dot(u)\n", + " return loss + pen\n", + "\n", + " # define initial guess\n", + " x0 = np.ones((T,)) #np.zeros((T,))\n", + "\n", + " # define optimization method\n", + " # By default, it will be 'BFGS', 'L-BFGS-B', or 'SLSQP' depending on the constraints and bounds\n", + " # If constraints, it will be 'SLSQP' (Sequential Least SQuares Programming)\n", + " # 'Nelder-Mead', 'Powell', 'CG', 'BFGS', 'Newton-CG', 'L-BFGS-B', 'TNC', 'COBYLA', 'SLSQP', 'dogleg', 'trust-ncg'\n", + " # 'Nelder-Mead', 'Powell', 'CG', 'Newton-CG', 'TNC', 'COBYLA', 'dogleg', 'trust-ncg' can not handle (eq) constraints\n", + " # 'BFGS', 'L-BFGS-B' do not work\n", + " method = 'SLSQP'\n", + "\n", + " # define 1st constraints: norm of 1\n", + " constraints = [{'type': 'eq', 'fun': lambda u: u.T.dot(u) - 1}]\n", + "\n", + " # optimize recursively\n", + " evals, evecs = [], []\n", + " messages = {}\n", + " for i in range(T):\n", + " if i != 0:\n", + " # add orthogonality constraint\n", + " constraints.append({'type': 'eq', 'fun': lambda u: u1.T.dot(u)})\n", + "\n", + " # minimize --> it returns an instance of OptimizeResult\n", + " u1 = minimize(f, x0, method=method, constraints=constraints)\n", + " if not u1.success:\n", + " messages[i] = u1.message\n", + "\n", + " # get 'eigenvalue'\n", + " evals.append(-u1.fun)\n", + "\n", + " # get solution\n", + " u1 = u1.x\n", + " evecs.append(u1)\n", + "\n", + " return np.array(evals), np.array(evecs).T, messages" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "evals, evecs, messages = pca_scipy(X)\n", + "P = evecs\n", + "print(messages)\n", + "print(evals)\n", + "print(P)\n", + "\n", + "# Plot\n", + "plt.figure(figsize=(10,4))\n", + "plt.subplot(1,2,1)\n", + "plt.title('PCA: eigenvalues')\n", + "plt.bar(np.array([0.,0.1]), evals, width=0.1)\n", + "plt.xlim(0.,1.)\n", + "\n", + "plt.subplot(1,2,2)\n", + "plt.title('PCA: data and eigenvectors')\n", + "plt.scatter(X[:,0], X[:,1], color='b')\n", + "plt.arrow(0, 0, P[0,0], P[1,0], length_includes_head = True, head_width = 0.15, color='b')\n", + "plt.arrow(0, 0, P[0,1], P[1,1], length_includes_head = True, head_width = 0.15, color='g')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Using CVXPY\n", + "\n", + "Note: You **cannot** use `cvxpy` for this problem, as we are trying to maximize a convex function, and `cvxpy` only accepts to maximize a concave fct." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import cvxpy as cvx\n", + "\n", + "# cache the 'covariance' matrix\n", + "C = X.T.dot(X)/(N-1)\n", + "\n", + "# define vector to optimize\n", + "u1 = cvx.Variable(X.shape[1])\n", + "print(cvx.quad_form(u1, C).is_dcp())\n", + "print(cvx.quad_form(u1, C).is_quadratic())\n", + "\n", + "# define objective fct to maximize\n", + "#f = cvx.Maximize(u1.T*C*u1) \n", + "f = cvx.Maximize(cvx.quad_form(u1, C)) # this does not work!\n", + "#f = cvx.Minimize(cvx.quad_form(u1, C)) # this works (if no constraints) but that is not what we want to achieve!\n", + "constraints = [u1.T*u1 == 1]\n", + "prob = cvx.Problem(f, constraints)\n", + "\n", + "result = prob.solve()\n", + "print(u1.value)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Using NLopt\n", + "\n", + "Nonlinear optimization algorithms that can handle nonlinear inequality and **equality** constraints are:\n", + "* ISRES (Improved Stochastic Ranking Evolution Strategy) $\\rightarrow$ global derivative-free\n", + "* COBYLA (Constrained Optimization BY Linear Approximations) $\\rightarrow$ local derivative-free\n", + "* SLSQP (Sequential Least-SQuares Programming) $\\rightarrow$ local gradient-based\n", + "* AUGLAG (AUGmented LAGrangian) $\\rightarrow$ global/local derivative-free/gradient based (determined based on the subsidiary algo)\n", + "\n", + "More information about the various algorithms can be found on this [link](https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Playground with NLopt\n", + "import nlopt\n", + "\n", + "nlopt_results = {1: 'success', 2: 'stop_val reached', 3: 'ftol reached', 4: 'xtol reached',\n", + " 5: 'maxeval reached', 6: 'maxtime reached', -1: 'failure', -2: 'invalid args',\n", + " -3: 'out of memory', -4: 'roundoff limited', -5: 'forced stop'}\n", + "n_iter = 0\n", + "N = X.shape[0]\n", + "\n", + "# cache the 'covariance' matrix\n", + "C = X.T.dot(X)/(N-1)\n", + "\n", + "# define random seed\n", + "nlopt.srand(125)\n", + "\n", + "# define which solver to use\n", + "#optimizer = nlopt.GN_ISRES\n", + "#optimizer = nlopt.LN_COBYLA\n", + "optimizer = nlopt.LD_SLSQP\n", + "#optimizer = nlopt.LD_AUGLAG # nlopt.AUGLAG, nlopt.AUGLAG_EQ, nlopt.LD_AUGLAG,\n", + " # nlopt.LD_AUGLAG_EQ, nlopt.LN_AUGLAG, nlopt.LN_AUGLAG_EQ\n", + "\n", + "# if nlopt.AUGLAG, we have to define a subsidiary algo\n", + "suboptimizer = nlopt.LD_SLSQP #nlopt.LN_COBYLA\n", + "\n", + "# define objective function\n", + "def f(x, grad):\n", + " global n_iter\n", + " n_iter += 1\n", + " if grad.size > 0:\n", + " grad[:] = 2*x.T.dot(C)\n", + " return x.T.dot(C).dot(x)\n", + "\n", + "# define norm constraint\n", + "def c1(x, grad):\n", + " if grad.size > 0:\n", + " grad[:] = 2*x\n", + " return (x.T.dot(x) - 1)\n", + "\n", + "# create optimizer\n", + "n = X.shape[1] # nb of parameters to optimize, size of the vector\n", + "opt = nlopt.opt(optimizer, n)\n", + "print(\"Algo: %s\" % opt.get_algorithm_name())\n", + "opt.set_max_objective(f)\n", + "\n", + "# if nlopt.GN_ISRES, we can define the population size\n", + "opt.set_population(0) # by default for ISRES: pop=20*(n+1)\n", + "\n", + "# if nlopt.AUGLAG, set the subsidiary algo\n", + "subopt = nlopt.opt(suboptimizer, n)\n", + "subopt.set_lower_bounds(-1)\n", + "subopt.set_upper_bounds(1)\n", + "#subopt.set_ftol_rel(1e-2)\n", + "#subopt.set_maxeval(100)\n", + "opt.set_local_optimizer(subopt)\n", + "\n", + "# define bound constraints (should be between -1 and 1 because the norm should be 1)\n", + "opt.set_lower_bounds(-1.)\n", + "opt.set_upper_bounds(1.)\n", + "\n", + "# define equality constraints\n", + "opt.add_equality_constraint(c1, 0)\n", + "#opt.add_equality_mconstraint(constraints, tol)\n", + "\n", + "# define stopping criteria\n", + "#opt.set_stopval(stopval)\n", + "opt.set_ftol_rel(1e-8)\n", + "#opt.set_xtol_rel(1e-4)\n", + "opt.set_maxeval(100000) # nb of iteration\n", + "opt.set_maxtime(10) # time in secs\n", + "\n", + "# define initial value\n", + "#x0 = np.zeros((n,))\n", + "x0 = np.array([0.1,0.1])\n", + "\n", + "# optimize\n", + "x = x0\n", + "try:\n", + " x = opt.optimize(x0)\n", + "except nlopt.RoundoffLimited as e:\n", + " pass\n", + "max_value = opt.last_optimum_value()\n", + "result = opt.last_optimize_result()\n", + "\n", + "print(\"Max value: %f\" % max_value)\n", + "print(\"Nb of iterations: %d\" % n_iter)\n", + "print(\"Result: %s\" % nlopt_results[result])\n", + "print(\"Optimized array:\")\n", + "print(x)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define PCA optimization method formally using NLopt\n", + "def center_data(X, normalize=False, copy=True):\n", + " if copy:\n", + " X = np.copy(X)\n", + "\n", + " # center the data\n", + " mean = X.mean(axis=0)\n", + " X -= mean\n", + " N = X.shape[0]\n", + " T = X.shape[1]\n", + "\n", + " # normalize\n", + " if normalize:\n", + " X /= X.std(axis=0)\n", + " \n", + " return X\n", + "\n", + "class OrthogonalConstraint(object):\n", + " \n", + " def __init__(self, v):\n", + " self.v = np.copy(v)\n", + " \n", + " def constraint(self, x, grad):\n", + " if grad.size > 0:\n", + " grad[:] = self.v\n", + " return (x.T.dot(self.v))\n", + " \n", + "\n", + "def pca_nlopt(X, method=None, submethod=None, rough_param=0.0, normalize=False, copy=True):\n", + " \"\"\"\n", + " Compute PCA on the given data using nlopt.\n", + " \n", + " :param (str) method: it can take the following value: 'SLSQP', 'ISRES',\n", + " 'COBYLA', 'AUGLAG'. By default, it will be 'SLSQP'.\n", + " :param (str) submethod: this needs to be defined only if method is 'AUGLAG'.\n", + " By default, it will be 'SLSQP'.\n", + " \"\"\"\n", + " # center the data\n", + " X = center_data(X, normalize=normalize, copy=copy)\n", + "\n", + " # compute 'covariance' matrix and cache it\n", + " N = X.shape[0]\n", + " C = X.T.dot(X) / (N-1)\n", + " \n", + " # define useful variables\n", + " nlopt_results = {1: 'success', 2: 'stop_val reached', 3: 'ftol reached', 4: 'xtol reached',\n", + " 5: 'maxeval reached', 6: 'maxtime reached', -1: 'failure', -2: 'invalid args',\n", + " -3: 'out of memory', -4: 'roundoff limited', -5: 'forced stop'}\n", + " n = X.shape[1] # nb of parameters to optimize, size of the vector\n", + " \n", + " # define random seed\n", + " nlopt.srand(125)\n", + "\n", + " # define which solver (and possibly subsolver) to use\n", + " def get_opt(method):\n", + " if method == 'ISRES':\n", + " return nlopt.opt(nlopt.GN_ISRES, n)\n", + " elif method == 'COBYLA':\n", + " return nlopt.opt(nlopt.LN_COBYLA, n)\n", + " elif method == 'SLSQP':\n", + " return nlopt.opt(nlopt.LD_SLSQP, n)\n", + " elif method == 'AUGLAG':\n", + " return nlopt.opt(nlopt.AUGLAG, n)\n", + " else:\n", + " raise NotImplementedError(\"The given method has not been implemented\")\n", + "\n", + " if method is None:\n", + " method = 'SLSQP' \n", + " opt = get_opt(method)\n", + " if method == 'AUGLAG':\n", + " if submethod is None:\n", + " submethod = 'SLSQP'\n", + " elif submethod == 'AUGLAG':\n", + " raise ValueError(\"Submethod should be different from AUGLAG\")\n", + " subopt = get_opt(submethod)\n", + " subopt.set_lower_bounds(-1)\n", + " subopt.set_upper_bounds(1)\n", + " #subopt.set_ftol_rel(1e-2)\n", + " #subopt.set_maxeval(100)\n", + " opt.set_local_optimizer(subopt)\n", + " \n", + " # define objective function\n", + " def f(x, grad):\n", + " if grad.size > 0:\n", + " grad[:] = 2*x.T.dot(C)\n", + " return x.T.dot(C).dot(x)\n", + "\n", + " # define norm constraint\n", + " def c1(x, grad):\n", + " if grad.size > 0:\n", + " grad[:] = 2*x\n", + " return (x.T.dot(x) - 1)\n", + " \n", + " # define objective function\n", + " opt.set_max_objective(f)\n", + " \n", + " # if nlopt.GN_ISRES, we can define the population size\n", + " opt.set_population(0) # by default for ISRES: pop=20*(n+1)\n", + " \n", + " # define bound constraints (should be between -1 and 1 because the norm should be 1)\n", + " opt.set_lower_bounds(-1.)\n", + " opt.set_upper_bounds(1.)\n", + "\n", + " # define equality constraints\n", + " opt.add_equality_constraint(c1, 0)\n", + " #opt.add_equality_mconstraint(constraints, tol)\n", + "\n", + " # define stopping criteria\n", + " #opt.set_stopval(stopval)\n", + " opt.set_ftol_rel(1e-8)\n", + " #opt.set_xtol_rel(1e-4)\n", + " opt.set_maxeval(100000) # nb of iteration\n", + " opt.set_maxtime(2) # time in secs\n", + "\n", + " # define initial value\n", + " x0 = np.array([0.1]*n) # important that the initial value ≠ 0 for the computation of the grad!\n", + "\n", + " evals, evecs, msgs = [], [], {}\n", + " for i in range(n):\n", + " if i > 0:\n", + " c = OrthogonalConstraint(x)\n", + " opt.add_equality_constraint(c.constraint, 0)\n", + " # optimize\n", + " try:\n", + " x = opt.optimize(x0)\n", + " except nlopt.RoundoffLimited as e:\n", + " pass\n", + "\n", + " # save values\n", + " evecs.append(x)\n", + " evals.append(opt.last_optimum_value())\n", + " msgs[i] = nlopt_results[opt.last_optimize_result()]\n", + "\n", + " return np.array(evals), np.array(evecs).T, msgs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "method = 'SLSQP' # 'SLSQP', 'COBYLA', 'ISRES', 'AUGLAG'\n", + "submethod = None\n", + "\n", + "evals, P, msgs = pca_nlopt(X, method=method, submethod=submethod)\n", + "print(msgs)\n", + "print(evals)\n", + "print(P)\n", + "\n", + "# Plot\n", + "plt.figure(figsize=(10,4))\n", + "plt.subplot(1,2,1)\n", + "plt.title('PCA: eigenvalues')\n", + "plt.bar(np.array([0.,0.1]), evals, width=0.1)\n", + "plt.xlim(0.,1.)\n", + "\n", + "plt.subplot(1,2,2)\n", + "plt.title('PCA: data and eigenvectors')\n", + "plt.scatter(X[:,0], X[:,1], color='b')\n", + "plt.arrow(0, 0, P[0,0], P[1,0], length_includes_head = True, head_width = 0.15, color='b')\n", + "plt.arrow(0, 0, P[0,1], P[1,1], length_includes_head = True, head_width = 0.15, color='g')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For comparison purpose, we obtained the following values with scipy ('SLSQP'):\n", + "\n", + "[ 7.93954312 0.06045688]
\n", + "[[ 0.83849224 -0.54491355]
\n", + " [ 0.54491353 0.83849226]]\n", + "\n", + "And these values using std PCA:\n", + "\n", + "[ 7.93954312 0.06045688]
\n", + "[[-0.83849224 0.54491354]
\n", + " [-0.54491354 -0.83849224]]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Using IPopt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Playground with IPopt\n", + "import ipopt\n", + "\n", + "# define useful vars\n", + "n = X.shape[1]\n", + "N = X.shape[0]\n", + "C = X.T.dot(X) / (N-1)\n", + "\n", + "# define initial value\n", + "x0 = np.array([0.1]*n)\n", + "\n", + "# define (lower and upper) bound constraints\n", + "lb = [-1]*n\n", + "ub = [1]*n\n", + "\n", + "# define constraints\n", + "cl = [1]\n", + "cu = [1]\n", + "\n", + "class Opt(object):\n", + " \n", + " def __init__(self, verbose=True):\n", + " self.verbose = verbose\n", + " self.iter_count = 0\n", + "\n", + " def objective(self, x):\n", + " # objective fct to minimize\n", + " return -x.T.dot(C).dot(x)\n", + " \n", + " def gradient(self, x):\n", + " # grad of the objective fct\n", + " return -2*x.T.dot(C)\n", + " \n", + " def constraints(self, x):\n", + " # norm constraint\n", + " return x.T.dot(x)\n", + " \n", + " def jacobian(self, x):\n", + " return 2*x\n", + " \n", + " #def hessian(self, x):\n", + " # pass\n", + " \n", + " def intermediate(self, alg_mod, iter_count, obj_value, inf_pr, inf_du, mu, d_norm,\n", + " regularization_size, alpha_du, alpha_pr, ls_trials):\n", + " if self.verbose:\n", + " print(\"Objective value at iteration #%d: %g\" % (iter_count, obj_value))\n", + " self.iter_count = iter_count\n", + "\n", + "opt = Opt(verbose=False)\n", + "nlp = ipopt.problem(n=n, m=len(cl), problem_obj=opt, lb=lb, ub=ub, cl=cl, cu=cu)\n", + "\n", + "x, info = nlp.solve(x0)\n", + "print(\"Max value: %f\" % -info['obj_val'])\n", + "print(\"Nb of iterations: %d\" % opt.iter_count)\n", + "print(\"Result: %s\" % info['status_msg'])\n", + "print(\"Optimized array:\")\n", + "print(x)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define PCA optimization method formally using IPopt\n", + "class NormConstraint(object):\n", + " \n", + " def __init__(self):\n", + " pass\n", + " \n", + " def constraint(self, x):\n", + " return x.T.dot(x)\n", + " \n", + " def jacobian(self, x):\n", + " return 2*x\n", + " \n", + "class OrthogonalConstraint(object):\n", + " \n", + " def __init__(self, v):\n", + " self.v = np.copy(v)\n", + " \n", + " def constraint(self, x):\n", + " return x.T.dot(self.v)\n", + " \n", + " def jacobian(self, x):\n", + " return self.v\n", + "\n", + "class IPopt(object):\n", + "\n", + " def __init__(self, verbose=True):\n", + " self.verbose = verbose\n", + " self.iter_count = 0\n", + " self.cons = []\n", + " \n", + " def add_constraint(self, c):\n", + " self.cons.append(c)\n", + "\n", + " def objective(self, x):\n", + " # objective fct to minimize\n", + " return -x.T.dot(C).dot(x)\n", + " \n", + " def gradient(self, x):\n", + " # grad of the objective fct\n", + " return -2*x.T.dot(C)\n", + " \n", + " def constraints(self, x):\n", + " return np.array([c.constraint(x) for c in self.cons])\n", + " \n", + " def jacobian(self, x):\n", + " return np.array([c.jacobian(x) for c in self.cons])\n", + " \n", + " #def hessian(self, x):\n", + " # pass\n", + " \n", + " def intermediate(self, alg_mod, iter_count, obj_value, inf_pr, inf_du, mu, d_norm,\n", + " regularization_size, alpha_du, alpha_pr, ls_trials):\n", + " if self.verbose:\n", + " print(\"Objective value at iteration #%d: %g\" % (iter_count, obj_value))\n", + " self.iter_count = iter_count\n", + "\n", + "def pca_ipopt(X, rough_param=0.0, normalize=False, copy=True):\n", + " \"\"\"\n", + " Compute PCA on the given data using ipopt.\n", + " \"\"\"\n", + " # center the data\n", + " X = center_data(X, normalize=normalize, copy=copy)\n", + " \n", + " # define useful vars\n", + " n = X.shape[1]\n", + " N = X.shape[0]\n", + " C = X.T.dot(X) / (N-1)\n", + "\n", + " # define initial value\n", + " x0 = np.array([0.1]*n)\n", + "\n", + " # define (lower and upper) bound constraints\n", + " lb = [-1]*n\n", + " ub = [1]*n\n", + "\n", + " # define constraints\n", + " cl = [1] + [0]*(n-1)\n", + " cu = [1] + [0]*(n-1)\n", + "\n", + " evals, evecs, msgs = [], [], {}\n", + " opt = IPopt(verbose=False)\n", + " opt.add_constraint(NormConstraint())\n", + " for i in range(n):\n", + " if i > 0:\n", + " opt.add_constraint(OrthogonalConstraint(x))\n", + " \n", + " i1 = i+1\n", + " nlp = ipopt.problem(n=n, m=len(cl[:i1]), problem_obj=opt,\n", + " lb=lb, ub=ub, cl=cl[:i1], cu=cu[:i1])\n", + " \n", + " # solve problem\n", + " x, info = nlp.solve(x0)\n", + " \n", + " evecs.append(x)\n", + " evals.append(-info['obj_val'])\n", + " msgs[i] = info['status_msg']\n", + " \n", + " return np.array(evals), np.array(evecs).T, msgs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "evals, P, msgs = pca_ipopt(X)\n", + "print(msgs)\n", + "print(evals)\n", + "print(P)\n", + "\n", + "# Plot\n", + "plt.figure(figsize=(10,4))\n", + "plt.subplot(1,2,1)\n", + "plt.title('PCA: eigenvalues')\n", + "plt.bar(np.array([0.,0.1]), evals, width=0.1)\n", + "plt.xlim(0.,1.)\n", + "\n", + "plt.subplot(1,2,2)\n", + "plt.title('PCA: data and eigenvectors')\n", + "plt.scatter(X[:,0], X[:,1], color='b')\n", + "plt.arrow(0, 0, P[0,0], P[1,0], length_includes_head = True, head_width = 0.15, color='b')\n", + "plt.arrow(0, 0, P[0,1], P[1,1], length_includes_head = True, head_width = 0.15, color='g')\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.12" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/tutorials/math/README.md b/tutorials/math/README.md new file mode 100644 index 0000000..8c5c32d --- /dev/null +++ b/tutorials/math/README.md @@ -0,0 +1,14 @@ +## Tutorials about mathematics + +This folder contains tutorials about mathematics used in machine learning or robotics. It might happen that sometimes I will redirect to other webpages that explain the concept better than me. + +* linear algebra +* multivariate calculus +* probability and statistics +* information theory +* variational inference +* topologies and manifolds +* tensor calculus +* differential geometry +* information geometry + diff --git a/tutorials/python/README.md b/tutorials/python/README.md new file mode 100644 index 0000000..9ed8d7e --- /dev/null +++ b/tutorials/python/README.md @@ -0,0 +1,65 @@ +## Python Tutorials + +Hey! So you want to learn Python? Good :thumbsup: you are at the right place! :smile: + +Here is what you have to know: + +- There are mainly 2 versions of Python: Python 2 and Python 3. The most known version of the former is Python 2.7, and is currently installed by default on Linux and MacOSX systems. You can try it by typing `python` in the terminal. To exit, push `Ctrl+D` or type `exit()`. Note that Python 2.7 will stop to be maintained in 2020. Now, I am pretty sure that you are asking yourself "what is the difference between Python 2.7 and Python 3?". Well, to be short "Python 2.x is legacy, Python 3.x is the present and future of the language" as mentioned [here](https://wiki.python.org/moin/Python2orPython3). However, if you are using ROS/Gazebo, I would still recommend to use Python 2.7 for now. Pratically, these 2 versions are similar except for few details, the `print` statements where you don't have to put parenthesis in Python 2.7, and for few specific libraries. Personally, most of my code works in Python 3 and Python 2.7. + +- Check first this ["Tutorial: Learn Python in 10 Minutes"](https://www.stavros.io/tutorials/python/) to get a general understanding of the language. It assumes that you already know at least one programming language such as C++, Java or Matlab. Otherwise, it will take you more than 10minutes :stuck_out_tongue_winking_eye: + +- Then check this [one](https://learnxinyminutes.com/docs/python/) for Python 2.7 or this [one](https://learnxinyminutes.com/docs/python3/) for Python 3. + +- Well now, you just have to practice! :grin: + +## About Libraries/Modules + +The libraries (aka modules) that are useful for roboticists and/or machine learning engineers are: +* `numpy` for vectors, matrices, and tensors. Here is the [tutorial](https://docs.scipy.org/doc/numpy-dev/user/quickstart.html). +* `matplotlib` for plotting. + * check [here](https://matplotlib.org/gallery.html) to see the gallery (basically it shows you what you can do with it). + * check [this tutorial](https://matplotlib.org/users/beginner.html) to learn on how to use it. The first link `Pyplot tutorial` should be enough to have a basic understanding. + * if you are interested by 3D plots, check this [one](https://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html) + * and if you want animations, look at this [blog (and the webpages it links to)](https://jakevdp.github.io/blog/2012/08/18/matplotlib-animation-tutorial/) +* `scipy` for scientific and mathematical tools such as optimization, interpolation, signal processing, etc. You can find the tutorials at the following [link](https://docs.scipy.org/doc/scipy/reference/tutorial/). + * If you are not satisfied with `scipy.optimize.minimize`, you can have a look at [`cvxpy`](http://www.cvxpy.org/en/latest/) (which uses different solvers such as [`cvxopt`](http://cvxopt.org/)) for convex optimization. For QP problems with cvxopt, you can have a look at this small [tutorial](https://courses.csail.mit.edu/6.867/wiki/images/a/a7/Qp-cvxopt.pdf). For general nonlinear optimization, you can check [`nlopt`](https://nlopt.readthedocs.io/en/latest/), or [`ipopt`](https://pypi.python.org/pypi/ipopt) with the corresponding [documentation](http://pythonhosted.org/ipopt/). Most of these softwares are written in C/C++/Fortran but provides python wrappers to call the various functions. +* `sklearn` for machine learning algorithms. You can have a look at this library with the tutorials and documentations at the following [link](http://scikit-learn.org/stable/). +* `pandas` when working with tables and data structures. It also provides you data analysis tools. Check [here](https://pandas.pydata.org/pandas-docs/stable/10min.html) for a "10min" tutorial. + +More specific libraries: +* Deep Learning (or if you want to work with Tensors, or use automatic derivation tools) + * `TensorFlow` (aka `TF`) which can be found [here](https://www.tensorflow.org/). + * `PyTorch` at the following [webpage](http://pytorch.org/). + * Difference between `PyTorch` and `TF`? You can have a look at this [blog](https://awni.github.io/pytorch-tensorflow/) and this [medium post](https://medium.com/towards-data-science/pytorch-vs-tensorflow-spotting-the-difference-25c75777377b). Quickly, TF is being developed by Google, is more stable, has a bigger community, deals very nicely with the hardware part, has a syntax close to numpy, and is the current tool to use if you want to develop a software product. Meanwhile, PyTorch has a better integration with Python, allows you to use dynamic graphs (instead of static ones as in TF) and experiment new ideas faster. Hovewer, it is still in its early phase and the syntax used is currently different from numpy. + * `Keras`: this library is in the process of being integrated into TensorFlow. Quick description: it provides higher functionalities and is used on top of Theano or TensorFlow. + * `Theano`: this library is no longer maintained. +* Gaussian Processes + * `GPy`: the repo can be found [here](https://github.com/SheffieldML/GPy). Note that `sklearn` also allows you to use basic GPs but it is not as complete as `GPy`. For tutorials, have a look [here](http://nbviewer.jupyter.org/github/SheffieldML/notebook/blob/master/GPy/index.ipynb) +* Reinforcement Learning + * RL Environments + * `OpenAI-Gym`: the repo is [here](https://github.com/openai/gym), the documentation can be found [here](https://gym.openai.com/docs/), and the environments are [here](https://gym.openai.com/envs/). +* Dynamic Movement Primitives + * `pydmps`: the repo can be found [here](https://github.com/studywolf/pydmps) and the tutorials (which are very nice) are on this [blog](https://studywolf.wordpress.com/category/robotics/dynamic-movement-primitive/). +* `rospy`: this library allows you to use [ROS](http://www.ros.org/) (the Robot Operating System). Check [here](http://wiki.ros.org/rospy) and [here](http://wiki.ros.org/rospy_tutorials) for the tutorials. + +C/C++ libraries with Python support: +* `OpenCV`: library for computer vision. +* `KDL`: kinematics and dynamics library. +* `RBDL`: rigid body dynamics library. + +## Nice tools + +Here are nice tools that you should know: +* [`Jupyter Notebook`](http://jupyter.org/). As mentioned on their webpage, "The Jupyter Notebook is an open-source web application that allows you to create and share documents that contain live code, equations, visualizations and explanatory text". +* [`pip`](https://pip.pypa.io/en/stable/) which is tool that allows you to install Python packages/modules/libraries. For instance, to install `numpy` you would type in the terminal: `pip install numpy`. +* [`Anaconda`](https://www.anaconda.com/) which allows you to create virtual environments and install packages like `pip`. +* [`Cython`](http://cython.org/): this one is not really a tool, but allows you to wrap C/C++ classes and functions, and use them in Python. You could also use it to write `Cython` code which would be faster than pure Python code. Having said that, it takes quite a bit of time to learn this language, and effectively use it. + +## For Matlab Users + +So you want to move from Matlab? Good, you are definitely at the right place :wink: +Here are few webpages that could be useful for you: +* The following [link](http://www.pyzo.org/python_vs_matlab.html) explains the difference between Matlab and Python. I would also add that Python has a bigger community (because it is free), and is becoming **the** programming language to use for machine learning. +* Numpy for Matlab Users: [Link1](http://mathesaurus.sourceforge.net/matlab-numpy.html) and [Link2](https://docs.scipy.org/doc/numpy-dev/user/numpy-for-matlab-users.html). + +That's all folks! :clap: diff --git a/tutorials/robotics/README.md b/tutorials/robotics/README.md new file mode 100644 index 0000000..190e19a --- /dev/null +++ b/tutorials/robotics/README.md @@ -0,0 +1,27 @@ +## Robotics tutorials + +This folder contains tutorials about robotics in the form of python notebooks and webpages, and are based on books, slides, videos and notes [1, 2, 3, 4, 5]. I hope these tutorials will be useful for people learning about robotics, and a possible supplement to their theoretical courses. + + +## References + +Here are the main references that I used while writing the tutorials: +1. "Robotics: Modelling, Planning and Control", Siciliano et al., 2010 +2. [CS223A - Introduction to Robotics](https://see.stanford.edu/course/cs223a) by Prof. Khatib. +3. [Course materials](http://www.diag.uniroma1.it/deluca/Teaching.php) from Prof. Alessandro De Luca + - [Course material Robotics 1](http://www.diag.uniroma1.it/deluca/rob1_en/material_rob1_en.html) + - [Course material Robotics 2](http://www.diag.uniroma1.it/deluca/rob2_en/material_rob2_en.html) +4. "Springer Handbook of Robotics", Siciliano et al., 2008 + +Note that each python notebook provides more accurate references. + + +## Citation + +If these tutorials were useful to you, please cite the corresponding references that I used and the pyrobolearn framework. + + +## Improvements? + +If you have any criticisms, questions, advices, or requests, don't hesitate to open an issue. + diff --git a/tutorials/robotics/cartpole.ipynb b/tutorials/robotics/cartpole.ipynb new file mode 100644 index 0000000..c741bcb --- /dev/null +++ b/tutorials/robotics/cartpole.ipynb @@ -0,0 +1,39 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Cartpole (inverted pendulum on a cart)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "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.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/tutorials/robotics/dynamics.ipynb b/tutorials/robotics/dynamics.ipynb new file mode 100644 index 0000000..65f595b --- /dev/null +++ b/tutorials/robotics/dynamics.ipynb @@ -0,0 +1,39 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Dynamics" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "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.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/tutorials/robotics/kinematics.ipynb b/tutorials/robotics/kinematics.ipynb new file mode 100644 index 0000000..79b12e8 --- /dev/null +++ b/tutorials/robotics/kinematics.ipynb @@ -0,0 +1,39 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Kinematics" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "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.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/tutorials/robotics/mass-spring-damper.ipynb b/tutorials/robotics/mass-spring-damper.ipynb new file mode 100644 index 0000000..3004917 --- /dev/null +++ b/tutorials/robotics/mass-spring-damper.ipynb @@ -0,0 +1,215 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Mass-spring-damper\n", + "\n", + "In this tutorial, we will describe the mechanics and control of the one degree of freedom translational mass-spring-damper system subject to a control input force. We will first derive the dynamic equations by hand. Then, we will derive them using the `sympy.mechanics` python package.\n", + "\n", + "The system on which we will work is depicted below:\n", + "\n", + "![mass-spring-damper system](http://ctms.engin.umich.edu/CTMS/Content/Introduction/System/Modeling/figures/mass_spring_damper.png)\n", + "\n", + "Note that in what follows, we use the notation $u(t) = F$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Mechanics" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Deriving the dynamical equations by hand\n", + "\n", + "#### 1.1 By using Newton equations\n", + "\n", + "Using Newton's law, we have:\n", + "\n", + "\\begin{align}\n", + " m \\ddot{x}(t) &= \\sum F_{ext} \\\\\n", + " &= - b \\dot{x}(t) - k x(t) + u(t)\n", + "\\end{align}\n", + "\n", + "#### 1.2 By using the Lagrange Method\n", + "\n", + "Let's first derive the kinematic and potential energies.\n", + "\n", + "\\begin{equation}\n", + " T = \\frac{1}{2} m \\dot{x} \\\\\n", + " V = - \\int \\vec{F} . \\vec{dl} = - \\int (-kx \\vec{1_x}) . dx \\vec{1_x} = \\frac{k x^2}{2}\n", + "\\end{equation}\n", + "\n", + "The Lagrangian is then given by:\n", + "\\begin{equation}\n", + " \\mathcal{L} = T - V = \\frac{1}{2} m \\dot{x} - \\frac{k x^2}{2}\n", + "\\end{equation}\n", + "\n", + "Using the Lagrange's equations we can derive the dynamics of the system:\n", + "\n", + "\\begin{equation}\n", + " \\frac{d}{dt} \\frac{\\partial \\mathcal{L}}{\\partial \\dot{q}} - \\frac{\\partial \\mathcal{L}}{\\partial q} = Q\n", + "\\end{equation}\n", + "\n", + "where $q$ are the generalized coordinates (in this case $x$), and $Q$ represents the non-conservative forces (input force, dragging or friction forces, etc).\n", + "\n", + "* $\\frac{d}{dt} \\frac{\\partial \\mathcal{L}}{\\partial \\dot{x}} = \\frac{d}{dt} m \\dot{x}(t) = m \\ddot{x}(t) $\n", + "* $\\frac{\\partial \\mathcal{L}}{\\partial x} = - k x(t) $\n", + "* $Q = - b \\dot{x}(t) + u(t) $\n", + "\n", + "which when putting everything back together gives us:\n", + "\n", + "\\begin{equation}\n", + " m \\ddot{x}(t) + b \\dot{x}(t) + k x(t) = u(t)\n", + "\\end{equation}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Deriving the dynamical equations using sympy" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "import sympy\n", + "import sympy.physics.mechanics as mechanics\n", + "from sympy import init_printing\n", + "init_printing(use_latex='mathjax')\n", + "from sympy import pprint" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "⎡ 2 ⎤\n", + "⎢ d d ⎥\n", + "⎢b⋅──(q(t)) + 1.0⋅k⋅q(t) + m⋅───(q(t)) - u(t)⎥\n", + "⎢ dt 2 ⎥\n", + "⎣ dt ⎦\n" + ] + } + ], + "source": [ + "# define variables \n", + "q = mechanics.dynamicsymbols('q')\n", + "dq = mechanics.dynamicsymbols('q', 1)\n", + "u = mechanics.dynamicsymbols('u')\n", + "\n", + "# define constants\n", + "m, k, b = sympy.symbols('m k b')\n", + "\n", + "# define the inertial frame\n", + "N = mechanics.ReferenceFrame('N')\n", + "\n", + "# define a particle for the mass\n", + "P = mechanics.Point('P')\n", + "P.set_vel(N, dq * N.x) # go in the x direction\n", + "Pa = mechanics.Particle('Pa', P, m)\n", + "\n", + "# define the potential energy for the particle (the kinematic one is derived automatically)\n", + "Pa.potential_energy = k * q**2 / 2.0\n", + "\n", + "# define the Lagrangian and the non-conservative force applied on the point P\n", + "L = mechanics.Lagrangian(N, Pa)\n", + "force = [(P, -b * dq * N.x + u * N.x)]\n", + "\n", + "# Lagrange equations \n", + "lagrange = mechanics.LagrangesMethod(L, [q], forcelist = force, frame = N)\n", + "pprint(lagrange.form_lagranges_equations())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Laplace transform and transfer function" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Applying the Laplace transform on the dynamic equation:\n", + "\n", + "\\begin{equation}\n", + " m \\ddot{x}(t) + b \\dot{x}(t) + k x(t) = u(t) \\stackrel{L}{\\rightarrow} m s^2 X(s) + b s X(s) + k X(s) = U(s)\n", + "\\end{equation}\n", + "\n", + "The transfer equation is given by:\n", + "\n", + "\\begin{equation}\n", + " H(s) = \\frac{X(s)}{U(s)} = \\frac{1}{m s^2 + b s + k}\n", + "\\end{equation}\n", + "\n", + "By calculating the pole:\n", + "\n", + "\\begin{equation}\n", + " m s^2 + b s + k = 0 \\Leftrightarrow s = \\frac{-b}{2m} \\pm \\sqrt{\\left(\\frac{b}{2m}\\right)^2 - \\frac{k}{m}}\n", + "\\end{equation}\n", + "\n", + "Note that $b, k, m > 0$ because they represent real physical quantities." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### LTI system\n", + "\n", + "We can rewrite the above equation as a first-order system of equations. Let's first define the state vector $\\pmb{x} = \\left[ \\begin{array}{c} x(t) \\\\ \\dot{x}(t) \\end{array} \\right]$ and the control vector $\\pmb{u} = \\left[ \\begin{array}{c} u(t) \\end{array} \\right]$, then we can rewrite the above equation in the form $\\pmb{\\dot{x}} = \\pmb{Ax} + \\pmb{Bu}$, as below:\n", + "\n", + "\\begin{equation}\n", + " \\left[ \\begin{array}{c} \\dot{x}(t) \\\\ \\ddot{x}(t) \\end{array} \\right] = \\left[ \\begin{array}{cc} 0 & 1 \\\\ -\\frac{k}{m} & -\\frac{b}{m} \\end{array} \\right] \\left[ \\begin{array}{c} x(t) \\\\ \\dot{x}(t) \\end{array} \\right] + \\left[ \\begin{array}{c} 0 \\\\ \\frac{1}{m} \\end{array} \\right] \\left[ \\begin{array}{c} u(t) \\end{array} \\right]\n", + "\\end{equation}\n", + "\n", + "If there is no $u(t)$, i.e. $u(t) = 0 \\; \\forall t$, then we have $\\pmb{\\dot{x}} = \\pmb{Ax}$. The solution to this system of equation is $\\pmb{x}(t) = e^{\\pmb{A}t} \\pmb{x}_0$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.12" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}