A1_vectors_matrices_in_numpy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Crash course 1: Vectors and matrices in numpy and scipy</center>\n",
    "### <center>Alfred Galichon (NYU & ScPo) and Clément Montes (ScPo)</center>\n",
    "## <center>'math+econ+code' masterclass on optimal transport and economic applications</center>\n",
    "#### <center>With python code examples</center>\n",
    "© 2018-2022 by Alfred Galichon. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274 are acknowledged, as well as inputs from contributors listed [here](http://www.math-econ-code.org/theteam).\n",
    "\n",
    "**If you reuse material from this masterclass, please cite as:**<br>\n",
    "Alfred Galichon, 'math+econ+code' masterclass on optimal transport and economic applications, January 2022. https://github.com/math-econ-code/mec_optim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introducing NumPy\n",
    "\n",
    "* Unlike R or Matlab, Python has no built-in matrix algebra interface. Fortunately, the NumPy library provides powerful matrix capabilities, on par with R or Matlab. Here is a quick introduction to vectorization, operations on vectors and matrices, higher-dimensional arrays, Kronecker products and sparse matrices, etc. in NumPy.\n",
    "\n",
    "* This is *not* a tutorial on Python itself. They are plenty good ones available on the web.\n",
    "\n",
    "* First, we load numpy (with its widely used alias):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In NumPy, an `array` is built from a lists as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = np.array([1,2,3])\n",
    "print(u)\n",
    "v = np.array([3,2,5])\n",
    "print(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can then add arrays as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.array([1,2,3])+np.array([3,2,5]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the difference between the + operator when applied to numpy arrays vs. when applied to lists:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[1,2,3]+[3,2,5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the latter case, it returns list concatenation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To input matrices in NumPy, one simply inputs a list of rows, which are themselves represented as lists."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[11,12],[21,22],[31,32]])\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `shape` attribute of an array indicated the dimension of that array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vectorization and memory order\n",
    "\n",
    "* Matrices in all mathematical softwares are represented in a *vectorized* way as a sequence of numbers in the computers memory. This representation can involve either stacking the lines, or stacking the columns.\n",
    "\n",
    "* Different programming languages can use either of the two stacking conventions:\n",
    "    + Stacking the lines (Row-major order) is used by `C`, and is the default convention for Python (NumPy). A matrix $M$ is represented by varying the last index first, i.e. a $2\\times2$ matrix will be represented as $vec_C\\left(M\\right) = \\left(M_{11}, M_{12}, M_{21}, M_{22}\\right).$ \n",
    "    + Stacking the columns (Column-major order) is used by `Fortran`, `Matlab`, `R`, and most underlying core linear algebra libraries (like BLAS). A 2x2x2 3-dimensional array $A$ will be represented by varying the first index first, then the second, i.e. $vec_C\\left(A\\right) = \\left( A_{111}, A_{112}, A_{121}, A_{122}, A_{211}, A_{212}, A_{221}, A_{222} \\right)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The command `flatten()` provides the vectorized representation of a matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember, NumPy represents matrices by **varying the last index first**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to reshape the matrix `a`, one modifies its `shape` attribute. The following reshapes the matrix `a` into a row vector. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.shape = 1,6\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The previous output evidences the fact that Python uses the row-major order: rows are stacked one after the other. \n",
    "To reshape the vector into a column vector, do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.shape = 6,1\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equivalently, one could have set `A.shape=6,-1`, where Python would replace `-1` by the integer needed for the formula to make sense (in this case, `1`). \n",
    "Another way to reshape is to use the method `reshape,` which returns a duplicate of the object with the requested shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A1=np.array(range(6))\n",
    "A2 = A1.reshape(3,2)\n",
    "print(\"A1=\\n\", A1)\n",
    "print(\"A2=\\n\",A2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `NumPy` also supports the column-major order, but you have to specifically ask for it, by passing the optional argument `order='F'`, where 'F' stands for `Fortran`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A3 = np.array(range(6)).reshape(3,2, order='F')\n",
    "A3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiplication "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multiplication of arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several ways to multiply two arrays using NumPy. The most commonly used is the following."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.ones((2,2))\n",
    "B = 3*np.eye(2)\n",
    "A@B #@ is left associative. If you have A@B@C, it will compute (A@B)@C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `np.matmul(A,B)` would give the same result as well, but it is more difficult to read `np.matmul(A,np.matmul(B,C))` than `A@B@C`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multiplication by a scalar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "4*np.eye(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above assignation of B corresponds to the multiplication by a scalar. It is the simplest broadcasting allowed by numpy (which makes this library more powerful than just using lists -it is also much quicker-). More on broadcasting will arrive later in that Notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kronecker product\n",
    "\n",
    "A very important identity is\n",
    "\\begin{align*}\n",
    "vec_C\\left(AXB\\right) = \\left(  A\\otimes B^\\top\\right)  vec_C\\left(X\\right),\n",
    "\\end{align*}\n",
    "where $vec_C$ is the vectorization under the C (row-major) order, and where the Kronecker product $\\otimes$ is defined as follows for 2x2 matrices (with obvious generalization):\n",
    "\n",
    "\\begin{align*}\n",
    "A\\otimes B=\n",
    "\\begin{pmatrix}\n",
    "a_{11}B & a_{12}B\\\\\n",
    "a_{21}B & a_{22}B\n",
    "\\end{pmatrix}.\n",
    "\\end{align*}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.eye(2)\n",
    "\n",
    "AXB = np.kron(A, B)\n",
    "print(\"A=\",A)\n",
    "print(\"B=\",B)\n",
    "print(\"AXB=\",AXB)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Type broadcasting in NumPy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The term broadcasting describes how NumPy treats arrays with different shapes during arithmetic operations. \n",
    "\n",
    "Subject to certain constraints, the smaller array is “broadcasted” across the larger array so that they have compatible shapes. Broadcasting provides a means of vectorizing array operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = 10*np.array([[1],[2],[3]]) #Simplest broadcasting\n",
    "B =  np.array([1,2])\n",
    "print('A=\\n',A)\n",
    "print('B=\\n',B)\n",
    "print('A+B=\\n',A+B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The operation `A[:,np.newaxis]` creates a new dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = np.array([3,4,5])\n",
    "print(v)\n",
    "print(v[:,np.newaxis])\n",
    "print(v[np.newaxis,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Arrays of larger dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a_3d_array = np.array([[[1, 2], [3, 4]], [[5, 6], [7, 8]]])\n",
    "a_3d_array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Standard functions can also support arrays with more than 2 dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a_multiarray = np.zeros((2,3,3,3))\n",
    "print(a_multiarray, a_multiarray.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Searching for a maximum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Maximum between 2 arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compare two arrays (say $x$ and $y$) component by component, it is convenient to use `np.maximum`. It returns an array $z$ such that $ \\forall i: z[i] = \\max(x[i],y[i])$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.maximum(np.array([2, 3, 4]), np.array([1, 5, 2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can even broadcast."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.maximum(np.eye(2), [0.5, 2]) # broadcasting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Highest component within an array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`np.max` and `np.argmax` respectively find the maximum entry of a given array along a specified axis, and its index. `np.min` and `np.argmin` perform similar functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[0, 1,3], [0, 5,7]])\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.max(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.argmax(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.max(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.min(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.argmin(axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If `axis` is not specified, the maximum will be taken over all the entries of the matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.max(A) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: if your array contains a nan, you can use `np.nanmax` in order to ignore those values while searching for the highest component."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summing all elements of an array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a similar fashion as above, `np.sum` sums the elements of an array over a given axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.sum( axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.sum(axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If `axis` is not specified, the sum is done over all the entries of the matrix.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sparse matrices in Scipy\n",
    "\n",
    "Sparse matrices are available in the `sparse` module of the `scipy` library. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.sparse as spr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 1000\n",
    "\n",
    "print('size of sparse identity matrix of size '+str(n) +' in MB = ' + str(spr.identity(n).data.size  / (1024**2)))\n",
    "\n",
    "print('size of dense identity matrix of size '+str(n) +' in MB  = ' + str(spr.identity(n).todense().nbytes  / (1024**2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Working with sparse matrices requires less storage. It is explained by the fact that while a dense matrix needs to encode every coefficient on a byte, sparse matrices only store the non-null coefficients. It is really convenient to work with such objects when it comes to matrices with really high sizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spr.identity(1000).data.size  , spr.identity(1000).todense().nbytes "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating sparse matrices..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ... with standard forms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "I5 = spr.identity(5)\n",
    "I5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can convert your sparse matrix into a dense one in order to visualise it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I5.todense()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ... from a dense matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a dense matrix and make it sparse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import uniform module to create random numbers\n",
    "from scipy.stats import uniform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(seed=42)\n",
    "dense_matrix = uniform.rvs(size=16, loc = 0, scale=2) #List of 16 random draws between 0 and 2\n",
    "dense_matrix = np.reshape(dense_matrix, (4, 4))\n",
    "dense_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dense_matrix[dense_matrix < 1] = 0 #Arbitrar criterion\n",
    "dense_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sparse_matrix = spr.csr_matrix(dense_matrix)\n",
    "print(sparse_matrix) #It prints a tuple giving the row and columns of the non-null component and its value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ... from scratch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can create two arrays containing respectively the rows and the column of the non-null coefficients.\n",
    "A third array would give the value of the non-null coefficient. The result is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# row indices\n",
    "row_ind = np.array([0, 1, 1, 3, 4])\n",
    "# column indices\n",
    "col_ind = np.array([0, 2, 4, 3, 4])\n",
    "# coefficients\n",
    "data = np.array([1, 2, 3, 4, 5], dtype=float)\n",
    "\n",
    "mat_coo = spr.coo_matrix((data, (row_ind, col_ind)))\n",
    "print(mat_coo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Every common operation seen below works with sparse matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I5 = spr.identity(5)\n",
    "I5 + np.ones((5,5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I5 + np.diag([1.,2.,3.,4.,5.])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I5 @ np.diag([1.,2.,3.,4.,5.])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kron_product = spr.kron(I5 , 10 *np.array([[1,2],[3,4]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kron_product.todense()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "A = np.ones((1000,1000))\n",
    "I5 = 3*spr.identity(1000)\n",
    "B = I5.todense()\n",
    "\n",
    "t0 = time.time()\n",
    "A@B\n",
    "t1 = time.time()\n",
    "A@I5\n",
    "t2 = time.time()\n",
    "\n",
    "print(t1-t0, t2-t1)"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}