07.1_Algorithms.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Algorithmic Complexity\n",
    "\n",
    "Profiling (e.g. with `timeit`) doesn’t tell us much about how an algorithm will perform on a different computer since it is determined by the hardware features. To compare performance in a device-indpendent fashion, a formalism (a.k.a the \"Big-O\") is used that characterizes functions in terms of their rates of growth as a function of the size *n* of the input.\n",
    "\n",
    "An algorithm is compared to a given function $g(n)$ with a well defined scaling with *n*, e.g. $n^2$; if the ratio of the two is bounded, than that algorithm is ${\\cal O}(g(n))$. Note that:\n",
    "* Only the largest terms in the scaling of $g(n)$ is kept in the notation\n",
    "* two algorithms can have the same complexity and have very different performance; the same complexity only implies that the difference in performance is independent of *n*.\n",
    "\n",
    "### Comparing bubble sort ${\\cal O}(n^2)$ and merge sort ${\\cal O}(n\\log{n})$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f1(n, k):\n",
    "    return k*n*n\n",
    "\n",
    "def f2(n, k):\n",
    "    return k*n*np.log(n)\n",
    "\n",
    "n = np.arange(0.1, 20001)\n",
    "plt.plot(n, f1(n, 1), c='blue')\n",
    "plt.plot(n, f2(n, 1000), c='red')\n",
    "plt.xlabel('Size of input (n)', fontsize=16)\n",
    "plt.ylabel('Number of operations', fontsize=16)\n",
    "plt.legend(['$\\mathcal{O}(n^2)$', '$\\mathcal{O}(n \\log n)$'], loc='best', fontsize=20);\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See [here](https://wiki.python.org/moin/TimeComplexity) for the complexity of operations on standard Python data structures. Note for instance that searching a list is much more expensive than searching a dictionary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Space Complexity\n",
    "\n",
    "We can aslo use ${\\cal O}$ notation in the same way to measure the space complexity of an algorithm.  The notion of space complexity becomes important when you data volume is of the same magntude or larger than the memory you have available. In that case, an algorihtm with high space complexity may end up having to swap memory constantly, and will perform far worse than its time complexity would suggest.\n",
    "\n",
    "Just as you should have a good idea of how your algorithm will scale with increasing *n*, you should also be able to know how much memroy your data structures will require. For example, if you had an $n\\times p$ matrix of integers, an $n\\times p$ matrix of flaots, and an $n\\times p$ matrix of complex floats, how large can $n$ and $p$ be before you run out of RAM to store them?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Notice how much overhead Python objects have\n",
    "# A raw integer should be 64 bits or 8 bytes only\n",
    "\n",
    "import sys\n",
    "print (sys.getsizeof(1))\n",
    "print (sys.getsizeof(1234567890123456789012345678901234567890))\n",
    "print (sys.getsizeof(3.14))\n",
    "print (sys.getsizeof(3j))\n",
    "print (sys.getsizeof('a'))\n",
    "print (sys.getsizeof('hello world'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (np.ones((100,100), dtype='byte').nbytes)\n",
    "print (np.ones((100,100), dtype='i2').nbytes)\n",
    "print (np.ones((100,100), dtype='int').nbytes) # default is 64 bits or 8 bytes\n",
    "print (np.ones((100,100), dtype='f4').nbytes)\n",
    "print (np.ones((100,100), dtype='float').nbytes) # default is 64 bits or 8 bytes\n",
    "print (np.ones((100,100), dtype='complex').nbytes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scipy: high-level scientific computing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `scipy` package contains various toolboxes dedicated to common issues in scientific computing. Its different submodules correspond to different applications, such as interpolation, integration, optimization, image processing, statistics, special functions, etc.\n",
    "\n",
    "`scipy` can be compared to other standard scientific-computing libraries, such as the GSL (GNU Scientific Library for C and C++), or Matlab’s toolboxes. `scipy` is the core package for scientific routines in Python; it is meant to operate efficiently on numpy arrays, so that numpy and `scipy` work hand in hand.\n",
    "\n",
    "Before implementing a routine, it is worth checking if the desired data processing is not already implemented in `scipy`. As non-professional programmers, scientists often tend to re-invent the wheel, which leads to buggy, non-optimal, difficult-to-share and unmaintainable code. By contrast, `scipy`’s routines are optimized and tested, and should therefore be used when possible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear algebra with scipy and numpy\n",
    "\n",
    "The `scipy.linalg` module provides standard linear algebra operations, relying on an underlying efficient implementation (BLAS, LAPACK).\n",
    "\n",
    "We will review a few examples and applications. Sometimes numpy implements those methods too: if a given algorithm is present both in numpy and scipy, typically the latter is performant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import linalg as la\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# to limit the printout\n",
    "%precision 4\n",
    "np.set_printoptions(suppress=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# norm of a vector\n",
    "v = np.array([1,2])\n",
    "print (la.norm(v))\n",
    "\n",
    "# distance between two vectors\n",
    "w = np.array([1,1])\n",
    "print (la.norm(v-w))\n",
    "\n",
    "# inner products \n",
    "print (v.dot(w))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Elaborate example: covariance matrix as outer product\n",
    "\n",
    "The inner product is just matrix multiplication of a 1×n vector with an n×1 vector.\n",
    "The outer product of two vectors is instead just the opposite. It is given by:\n",
    "\n",
    "$$\n",
    "v\\otimes w=vw^t\n",
    "$$\n",
    "\n",
    "Note that $v$ and $w$ are column vectors. The result of the inner product is a scalar. The result of the outer product is a matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (np.outer(v,w))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n, p = 10, 4\n",
    "v = np.random.random((p,n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.cov(v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = v - v.mean(1)[:, np.newaxis]\n",
    "w.dot(w.T)/(n - 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Traces and determinants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 6\n",
    "M = np.random.randint(100,size=(n,n))\n",
    "print(M,'\\n')\n",
    "print ('determinant:',la.det(M),'\\n')\n",
    "print ('trace:',M.trace(),'\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matrix Decomposition\n",
    "\n",
    "Often data analysis problems boil down to solving linear systems. An example is the [Netflix Competition](https://en.wikipedia.org/wiki/Netflix_Prize), where a matrix of $400000\\times18000$ (ratings times movies) needed to be dealt with. \n",
    "\n",
    "Matrix decompositions are an important step in solving linear systems in a computationally efficient manner.\n",
    "\n",
    "### Lower-Upper factorization\n",
    "\n",
    "Let A be a square matrix. An LU factorization refers to the factorization of A, with proper row and/or column orderings or permutations, into two factors – a lower triangular matrix L and an upper triangular matrix U:\n",
    "\n",
    "$A=LU$\n",
    "\n",
    "when solving a system of linear equations, $Ax=b=LUx$, the solution is done in two logical steps:\n",
    "1. solve $Ly=b$ for $y$.\n",
    "2. solve $Ux=y$ for $x$.\n",
    "\n",
    "Often a permutation $P$ is needed (*partial pivoting*) to best reorder the raws of the original matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1,3,4],[2,1,3],[4,1,2]])\n",
    "print(A,'\\n')\n",
    "\n",
    "P, L, U = la.lu(A)\n",
    "print(np.dot(P.T, A),'\\n')\n",
    "print(np.dot(L, U),'\\n')\n",
    "print(P,'\\n')\n",
    "print(L,'\\n')\n",
    "print(U,'\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Eigendecomposition\n",
    "\n",
    "Given an $n\\times n$ matrix $A$, with $\\det{A}\\ne0$, then there exist n  linearly independent eigenvectors and $A$ may be decomposed in the following manner:\n",
    "\n",
    "$$ \n",
    "A=V\\Lambda V^{-1}\n",
    "$$\n",
    "\n",
    "where $\\Lambda$ is a diagonal matrix whose diagonal entries are the eigenvalues of $A$ and the columns of $V$ are the corresponding eigenvectors of $A$.\n",
    "\n",
    "Eigenvalues are roots of the *characteristic polynomial* of $A$:\n",
    "\n",
    "$$\n",
    "\\det{A-\\lambda I}=0\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[0,1,1],[2,1,0],[3,4,5]])\n",
    "\n",
    "l, V = la.eig(A)\n",
    "print(np.dot(V,np.dot(np.diag(l), la.inv(V))),'\\n')\n",
    "print(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Singular Value Decomposition\n",
    "\n",
    "Another important matrix decomposition is singular value decomposition or SVD. For any $m\\times n$ matrix $A$, we may write:\n",
    "\n",
    "$$A=UDV^T$$\n",
    "\n",
    "where $U$ is a orthogonal $m\\times m$ matrix, $D$ (spectrum) is a rectangular, diagonal $m\\times n$ matrix with diagonal entries $d_1,\\dots,d_m$ all non-negative, $V$ is an orthogonal $n\\times n$ matrix.\n",
    "\n",
    "The singular-value decomposition is a generalization of the eigendecomposition in the sense that it can be applied to any $m \\times n$ matrix whereas eigenvalue decomposition can only be applied to diagonalizable matrices. \n",
    "\n",
    "Given an SVD of $A$, as described above, the following holds:\n",
    "\n",
    "$$\n",
    "A^T A = VD^TU^T UDV^T = VD^TDV^T \n",
    "$$\n",
    "$$\n",
    "A A^T = UD^TV^T VDU^T = UD^TDU^T \n",
    "$$\n",
    "\n",
    "The right-hand sides of these relations describe the eigenvalue decompositions of the left-hand sides. Consequently:\n",
    "* The columns of V (right-singular vectors) are eigenvectors of $A^TA$.\n",
    "* The columns of U (left-singular vectors) are eigenvectors of $AA^T$.\n",
    "* The non-zero elements of D (non-zero singular values) are the square roots of the non-zero eigenvalues of $A^TA$ or $AA^T$.\n",
    "\n",
    "A geometrical representation of SVD is given by the following figure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "Image(\"Singular-Value-Decomposition.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m, n = 5, 4\n",
    "A = np.random.randn(m, n) + 1.j*np.random.randn(m, n)\n",
    "print (A,'\\n')\n",
    "\n",
    "U, spectrum, Vt = la.svd(A)\n",
    "\n",
    "print(\"shapes:\", U.shape,  spectrum.shape, Vt.shape)\n",
    "\n",
    "print (spectrum,'\\n')\n",
    "print (U,'\\n')\n",
    "print (Vt,'\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = np.zeros((m, n))\n",
    "for i in range(min(m, n)):\n",
    "    D[i, i] = spectrum[i]\n",
    "SVD = np.dot(U, np.dot(D, Vt))\n",
    "np.allclose(SVD, A)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly scipy provide already a \"solve\" method for the linear systems of the type:\n",
    "\n",
    "$$A x = b $$\n",
    "\n",
    "still, knowing a little bit what are the algorithms underneath comes handy sometimes, e.g. the solve method can be instructed about what kind of matrix $A$ is likely to be (symmetric, hermitian, positive definite, etc.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[3, 2, 0], [1, -1, 0], [0, 5, 1]])\n",
    "b = np.array([2, 4, -1])\n",
    "x = la.solve(A, b)\n",
    "print (x)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.dot(A, x) == b\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Principal Component Analysis \n",
    "\n",
    "Principal Components Analysis (PCA) aims at finding and ranking all the eigenvalues and eigenvectors of a given dataset's covariance matrix. This is useful because high-dimensional data (with $p$ features) may have nearly all their variation in a small number of dimensions $k$, i.e. in the subspace spanned by the eigenvectors of the covariance matrix that have the $k$ largest eigenvalues. If we project the original data into this subspace, we can have a dimension reduction (from $p$ to $k$) with hopefully little loss of information.\n",
    "\n",
    "Numerically, PCA can be done either by means of eigendecomposition on the covariance matrix or via SVD on the data matrix. Even though the latter is usually preferred, let's have a look first at the former.\n",
    "\n",
    "To start with let's recall the definition of the covariance matrix (of 2 variables):\n",
    "\n",
    "$$\n",
    "{\\rm Cov}(X,Y)=\\frac{\\sum_{i=1}^n (X_i-\\bar{X})(Y_i-\\bar{Y})}{n-1}\n",
    "$$\n",
    "\n",
    "with Cov$(X,X)$ the variance of the variable $X$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# construct a dataset with a skewed 2D distribution\n",
    "mu = [0,0]\n",
    "sigma = [[0.6,0.2],[0.2,0.2]]\n",
    "n = 1000\n",
    "X = np.random.multivariate_normal(mu, sigma, n).T\n",
    "\n",
    "plt.scatter(X[0,:], X[1,:], alpha=0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the covariance matrix\n",
    "np.cov(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now find the eigenvectors of the covariance matrix..\n",
    "l, V = np.linalg.eig(np.cov(X))\n",
    "#l, V = la.eig(np.cov(X))\n",
    "\n",
    "# First recall that V is an orthogonal matrix (and thus its transpose is also its inverse)\n",
    "V.dot(V.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# ..and draw them (rescaling their module)\n",
    "plt.scatter(X[0,:], X[1,:], alpha=0.2)\n",
    "scale_factor=3\n",
    "for li, vi in zip(l, V.T):\n",
    "    plt.plot([0, scale_factor*li*vi[0]], [0, scale_factor*li*vi[1]], 'r-', lw=2)\n",
    "plt.axis([-3,3,-3,3])\n",
    "plt.title('Eigenvectors of covariance matrix scaled by eigenvalue.');\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the case the features of the datasets have all zero mean, the covariance matrix is of the form:\n",
    "\n",
    "$$\n",
    "{\\rm Cov}(X)=\\frac{XX^T}{n-1}\n",
    "$$\n",
    "\n",
    "and thus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l0, V0 = np.linalg.eig(np.dot(X, X.T)/(n-1))\n",
    "\n",
    "plt.scatter(X[0,:], X[1,:], alpha=0.2)\n",
    "for li, vi in zip(l0, V0.T):\n",
    "    plt.plot([0, scale_factor*li*vi[0]], [0, scale_factor*li*vi[1]], 'r-', lw=2)\n",
    "plt.axis([-3,3,-3,3]);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use the eigenvectors and eigenvalues to rotate the data, i.e. take the eigenvectors as new basis vectors and redefine the data points w.r.t this new basis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rotate all the data points accordingly to the new base\n",
    "Xp = np.dot(V0.T, X)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# then plot the rotated dataset and its \"axes\"\n",
    "plt.scatter(Xp[0,:], Xp[1,:], alpha=0.2)\n",
    "for li, vi in zip(l0, np.diag([1]*2)):\n",
    "    plt.plot([0, scale_factor*li*vi[0]], [0, scale_factor*li*vi[1]], 'r-', lw=2)\n",
    "plt.axis([-3,3,-3,3]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, if we only use the first column of $xp$, we will have the projection of the data onto the first principal component, capturing the majority of the variance in the data with a single featrue that is a linear combination of the original features.\n",
    "\n",
    "We may need to transform the (reduced) data set to the original feature coordinates for interpreation. This is simply another linear transform (matrix multiplication)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xpp = np.dot(V0, Xp)\n",
    "plt.scatter(Xpp[0,:], Xpp[1,:], alpha=0.2)\n",
    "for li, vi in zip(l0, V0.T):\n",
    "    plt.plot([0, scale_factor*li*vi[0]], [0, scale_factor*li*vi[1]], 'r-', lw=2)\n",
    "plt.axis([-3,3,-3,3]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dimension reduction via PCA\n",
    "\n",
    "\n",
    "Given the spectral decomposition:\n",
    "\n",
    "$$ \n",
    "A=V\\Lambda V^{-1}\n",
    "$$\n",
    "\n",
    "with $\\Lambda$ of rank $p$. Reducing the dimensionality to $k<p$ simply means setting to zero all bu first $k$ diagonal values (ordered from the largest to the smaller in module; that is the default in numpy/scipy).\n",
    "\n",
    "In this way we catch the most relevant part of its variability (covariance)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l, V = np.linalg.eig(np.cov(X))\n",
    "Lambda=np.diag(l)\n",
    "print (Lambda)\n",
    "print (\"A.trace():\", np.cov(X).trace())\n",
    "print (\"Lambda.trace():\", Lambda.trace())\n",
    "\n",
    "print (Lambda[0,0]/Lambda.trace())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the trace is invariant under change of basis, the total variability is also unchaged by PCA. By keeping only the first $k$ principal components, we can still “explain” \n",
    "$\\sum_1^k \\lambda_i/\\sum_1^p \\lambda_i$ of the total variability. Sometimes, the degree of dimension reduction is specified as keeping enough principal components so that (say) 90% fo the total variability is exlained."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SVD for PCA\n",
    "\n",
    "We saw that SVD is a decomposition of the data matrix $X=UDV^T$ where U and V are orthogonal matrices and D is a diagnonal matrix.\n",
    "\n",
    "\n",
    "Compare with the eigendecomposition of a matrix $A=W\\Lambda W^{−1}$, we see that SVD gives us the eigendecomposition of the matrix $XX^T$, which as we have just seen, is basically a scaled version of the covariance for a data matrix with zero mean, with the eigenvectors given by $U$ and eigenvalues by $D^2$ (scaled by n−1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U, spectrum, Vt = np.linalg.svd(X)\n",
    "\n",
    "l_svd = spectrum**2/(n-1)\n",
    "V_svd = U\n",
    "plt.scatter(X[0,:], X[1,:], alpha=0.2)\n",
    "for li, vi in zip(l_svd, V_svd):\n",
    "    plt.plot([0, scale_factor*li*vi[0]], [0, scale_factor*li*vi[1]], 'r-', lw=2)\n",
    "plt.axis([-3,3,-3,3]);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (\"eigendecomposition:\",l)\n",
    "print (\"SVD:\",l_svd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (\"eigendecomposition:\",V)\n",
    "print (\"SVD:\",V_svd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}