path: root/PVCM/cama/en/ma1 np05 Notation Einstein.ipynb

{
 "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,
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "## Introduction to Einstein Notation\n",
    "\n",
    "This notation allows for expressing an extraordinary number of operations concisely. This section is not\n",
    "necessary for the following but it is an interesting intellectual exercise.\n",
    "\n",
    "The basic idea is to sum the terms of an equation when the same index appears twice and is not defined elsewhere. Thus:\n",
    "\n",
    "$$ A_{i,i} \\quad \\textrm{means} \\quad \\sum_{i=1}^N A_{i,i} \\; \\textrm{(the trace of the matrix)}$$\n",
    "\n",
    "If we look at the matrix multiplication $A\\, B$, for every index (i,j) of the result we have:\n",
    "\n",
    "$$ C_{i,j} = A_{i,k} \\, B_{k,j} \\quad \\textrm{i.e.} \\quad  C_{i,j} = \\sum_{k=1}^N A_{i,k} \\, B_{k,j} $$ \n",
    "\n",
    "The full name of the Einstein notation being the Einstein summation convention, the Numpy function is [`einsum`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html). Here is how it works for our first 2 examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trace 'ii' :  12 \n",
      "\n",
      "Multiplication matricielle A A 'ij,jk->ik' :\n",
      "[[ 15  18  21]\n",
      " [ 42  54  66]\n",
      " [ 69  90 111]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "A = np.arange(9).reshape(3,3)\n",
    "\n",
    "print(\"Trace 'ii' : \", np.einsum('ii', A), '\\n')   # 0 + 4 + 8 = 12\n",
    "\n",
    "print(\"Multiplication matricielle A A 'ij,jk->ik' :\")\n",
    "print(np.einsum('ij,jk->ik', A, A))  # notez que j'ai nommé différement les indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 15,  18,  21],\n",
       "       [ 42,  54,  66],\n",
       "       [ 69,  90, 111]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.dot(A)  # on vérifie"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "It is noted that the arguments of `einsum` are:\n",
    "\n",
    "* the summation rule in a string with a comma to separate each component\n",
    "* the components on which the rule applies\n",
    "\n",
    "We can go a little further with Numpy. Here are all the rules that `einsum` uses:\n",
    "\n",
    "#### Basic and additional rules\n",
    "\n",
    "1. a repeated index implies summation over that index unless that index is mentioned in the result<br/>(see example of the diagonal of A below for the exception)\n",
    "2. indices that repeat from one component to another imply that the referenced elements will be multiplied together<br/> (see example of the matrix product)\n",
    "3. a letter omitted in the result (after `->`) implies summation over that index<br/> (see example of summing the elements of a vector below)\n",
    "4. if you don't put the arrow, einsum will put it for you with on the right all the indices that are not doubled arranged in alphabetical order<br/> (see example of the transpose below)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "Here is a list of operations taken from the [blog of Dr. Goulu](https://www.drgoulu.com/2016/01/17/einsum):\n",
    "\n",
    "<table align=\"center\">\n",
    "<tbody>\n",
    "<tr>\n",
    "<th>einsum Signature</th>\n",
    "<th>Numpy Equivalent</th>\n",
    "<th>Description</th>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('i-&gt;', v)</code></td>\n",
    "<td><code>sum(v)</code></td>\n",
    "<td>sum of values of vector v</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('i,i-&gt;i', u, v)</code></td>\n",
    "<td><code>u \\* v</code></td>\n",
    "<td>element-wise multiplication of vectors u and v</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('i,i', u, v)</code></td>\n",
    "<td><code>inner(u, v)</code></td>\n",
    "<td>dot product of u and v</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('i,j', u, v)</code></td>\n",
    "<td><code>outer(u, v)</code></td>\n",
    "<td>dyadic product of u and v</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ij', A)</code></td>\n",
    "<td><code>A</code></td>\n",
    "<td>returns matrix A</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ji', A)</code></td>\n",
    "<td><code>A.T</code></td>\n",
    "<td>transpose of A</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ii-&gt;i', A)</code></td>\n",
    "<td><code>diag(A)</code></td>\n",
    "<td>diagonal of A</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ii', A)</code></td>\n",
    "<td><code>trace(A)</code></td>\n",
    "<td>sum of the diagonal of A</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td class=\"tg-031e\"><code>('ij-&gt;', A)</code></td>\n",
    "<td class=\"tg-031e\"><code>sum(A)</code></td>\n",
    "<td class=\"tg-031e\">sum of values of A</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ij-&gt;j', A)</code></td>\n",
    "<td><code>sum(A, axis=0)</code></td>\n",
    "<td>sum of columns of A</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ij-&gt;i', A)</code></td>\n",
    "<td><code>sum(A, axis=1)</code></td>\n",
    "<td>sum of rows of A</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ij,ij-&gt;ij', A, B)</code></td>\n",
    "<td><code>A \\* B</code></td>\n",
    "<td>element-wise transposed matrix multiplication of A and B</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ij,jk', A, B)</code></td>\n",
    "<td><code>dot(A, B)</code></td>\n",
    "<td>dot product of A and B</td> </tr>\n",
    "<tr>\n",
    "<td><code>('ij,jk-&gt;ij', A, B)</code></td>\n",
    "<td><code>inner(A, B)</code></td>\n",
    "<td>inner product of A and B</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ij,jk-&gt;ijk', A, B)</code></td>\n",
    "<td><code>A[:, None] \\* B</code></td>\n",
    "<td>each row of A multiplied by B</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td><code>('ij,kl-&gt;ijkl', A, B)</code></td>\n",
    "<td><code>A[:, :, None, None] \\* B</code></td>\n",
    "<td>each value of A multiplied by B</td>\n",
    "</tr>\n",
    "</tbody>\n",
    "</table>\n",
    "\n",
    "The `None` in the last two lines is a way to reshape an array. So `np.arange(6)` is a 1-dimensional array, `np.arange(6)[:]` is the same 1-dimensional array, while `np.arange(6)[:,None]` is a 2-dimensional array namely `6 x 1` when `np.arange(6)[None,:,None]` has 3 dimensions: `1 x 6 x 1`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "## Practical Application\n",
    "\n",
    "We will compare the performance of `einsum` and the corresponding Numpy functions. To do this, calculate\n",
    "\n",
    "* the cube of each element of a vector\n",
    "* of a square matrix, $A^3$,\n",
    "\n",
    "with `einsum` and without. Compare the execution speed in all cases with 10000 elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a vous de jouer\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solution :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11.5 µs ± 112 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
      "15.2 µs ± 131 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
     ]
    }
   ],
   "source": [
    "u = np.random.random(10000)\n",
    "\n",
    "%timeit u*u*u\n",
    "%timeit np.einsum('i,i,i->i', u,u,u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "138 µs ± 9.48 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "134 ms ± 1.12 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "A = u.reshape(100,100)\n",
    "\n",
    "%timeit A.dot(A).dot(A)\n",
    "%timeit np.einsum('ij,jk,kl->il', A, A, A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "We observe that `einsum` is slower (\\*). But looking on the web, we see that this has not always been the case and that it is related to the version of the Numpy library. In conclusion, if we want performance, we should test our code beforehand to choose the fastest method.\n",
    "\n",
    "(\\*) slightly slower for vector calculation but 1000 times slower on my laptop for `A.dot(A)` (all my processors are running at 100% while the calculation by `einsum` is only done on 1 processor). The `A.dot(A)` version is much faster thanks to the MKL library used by Numpy on my machine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ]
}