{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting Started with NumPy and SciPy\n",
    "\n",
    "© Prof. Branislav K. Nikolić, University of Delaware\n",
    "\n",
    "[PHYS824: Nanophysics & Nanotechnology](https://wiki.physics.udel.edu/phys824) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[NumPy](https://numpy.org/) stands for Numerical Python. NumPy is a Python library used for working with arrays. It also has functions for working in domain of linear algebra, Fourier transform, and matrices. Although Python contains lists that serve the purpose of arrays, they are slow to process. NumPy aims to provide an array object that is up to 50x faster that traditional Python lists. The array object in NumPy is called ndarray, it provides a lot of supporting functions that make working with ndarray very easy.\n",
    "NumPy arrays are stored at one continuous place in memory unlike lists, so processes can access and manipulate them very efficiently. This behavior is called locality of reference in computer science. This is the main reason why NumPy is faster than lists. Also it is optimized to work with latest CPU architectures.\n",
    "\n",
    "[SciPy](https://www.scipy.org/) builds on NumPy and it contains high-level numerical routines dedicated to common issues in scientific computing. Its different submodules correspond to different applications, such as numerical integration, interpolation, optimization, Fast Fourier transform, signal and image processing, statistics, ordinary differential equation solvers, special functions, etc. 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. 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, routines within SciPy are optimized and tested, and should therefore be used when possible. The main Python package for linear algebra is SciPy subpackage [`scipy.linalg`](https://docs.scipy.org/doc/scipy/reference/linalg.html). \n",
    "\n",
    "Let us import both packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg as la"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For plotting we import matplotlib.pyplot package with option `%matplotlib inline`, so that plots are displayed in the notebook and not in a new window:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. NumPy arrays\n",
    "\n",
    "Let's begin with a quick review of [NumPy arrays](../../scipy/numpy/). We can think of a 1D NumPy array as a list of numbers. We can think of a 2D NumPy array as a matrix. And we can think of a 3D array as a cube of numbers.\n",
    "When we select a row or column from a 2D NumPy array, the result is a 1D NumPy array (called a [slice](https://docs.scipy.org/doc/numpy/reference/arrays.indexing.html#basic-slicing-and-indexing)). This is different from MATLAB where when you select a column from a matrix it is returned as a column vector which is a 2D MATLAB matrix. It can get a bit confusing and so we need to keep track of the shape, size and dimension of our NumPy arrays."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Interactive help"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.array?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.lookfor('create array')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.con*?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Array attributes\n",
    "\n",
    "Create a 1D (one-dimensional) NumPy array and verify its dimensions, shape and size. This apparently looks like a row vector in Matlab or bra vector $\\langle a|$ in the Dirac notation in quantum mechanics. But to mimic those properly those two concepts, take a look at [Section 5.2](#5.2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1,3,-2,1])\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Verify the number of dimensions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.ndim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Verify the shape of the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The shape of an array is returned as a [Python tuple](https://docs.python.org/3/tutorial/datastructures.html#tuples-and-sequences). The output in the cell above is a tuple of length 4. And we verify the size of the array (ie. the total number of entries in the array):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create column vector in Matlab or ket vector $|a \\rangle$ in quantum mechanics as $N \\times 1$ matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.array([[1],\n",
    "[3],\n",
    "[-2],\n",
    "[1]])\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b.ndim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b.size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a 2D (two-dimensional) NumPy array (i.e., matrix):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = np.array([[1,2],[3,7],[-1,5]])\n",
    "print(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Verify the number of dimensions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M.ndim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Verify the shape of the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, verify the total number of entries in the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M.size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select a column from a 2D NumPy array and we get a 1D array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col = M[:,1] \n",
    "print(col)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Verify the dimensions, shape and size of the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Dimensions:', col.ndim)\n",
    "print('Shape:', col.shape)\n",
    "print('Size:', col.size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we select a row of column from a 2D NumPy array, the result is a 1D NumPy array. However, we may want to select a column as a 2D column vector. This requires us to use the [reshape](https://docs.scipy.org/doc/numpy/reference/generated/numpy.reshape.html) method. For example, create a 2D column vector from the 1D slice selected from the matrix $M$ above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "column = np.array([2,7,5]).reshape(3,1)\n",
    "print(column)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Verify the dimensions, shape and size of the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Dimensions:', column.ndim)\n",
    "print('Shape:', column.shape)\n",
    "print('Size:', column.size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, variables `col` and `column` are different types of objects even though they have the \"same\" data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Functions for creating arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(1, 9, 2) # start, end (exclusive), step\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.linspace(0, 1, 5,  endpoint=True) # start, end, num-points\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.ones((3, 3))\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.zeros((3, 3))\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = np.eye(3)\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = np.diag(np.array([1, 2, 3, 4]))\n",
    "print(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Random numbers (Mersenne Twister PRNG):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.random.rand(5)  # uniform in [0, 1]\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.random.randn(5) # Gaussian\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1234) # setting the random seed "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4 Matrix visualization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = np.random.rand(50, 50)\n",
    "plt.imshow(image, cmap=plt.cm.hot)    \n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.5 Basic data types\n",
    "\n",
    "You may have noticed that, in some instances, array elements are displayed with a trailing dot (e.g. 2. vs 2). This is due to a difference in the data-type used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1, 2, 3])\n",
    "a.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.array([1., 2., 3.])\n",
    "b.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can explicitly specify which data-type you want:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = np.array([1, 2, 3], dtype=float)\n",
    "c.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **default** data type is floating point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.ones((3, 3))\n",
    "a.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = np.array([1+2j, 3+4j, 5+6*1j]) # complex with IMAGINARY UNIT denoted by 1j\n",
    "d.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "e = np.array([True, False, False, True]) # Boolean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = np.array(['Bonjour', 'Hello', 'Hallo'])\n",
    "f.dtype     # <--- strings containing max. 7 letters "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "## 2. Matrix operations and functions\n",
    "\n",
    "### 2.1 Arithmetic operations\n",
    "\n",
    "Arithmetic [array operations](../scipy/numpy/#operations-and-functions) `+`, `-`, `/`, `*` and `**` are performed **elementwise** on NumPy arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1,4],[9,16]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A * A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take square root of each element of matrix $A$ (which is not the same as $\\sqrt{A}$, see [Section 5.4](#5.4)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sqrt(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take exp of each element of matrix `np.sqrt(A)` (which is not the same as $\\exp{A}$, see [Section 6.4](#5.4)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.exp(np.sqrt(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = np.array([[3,4],[-1,5]])\n",
    "print(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But `np.sqrt` of $M$ gives nan:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sqrt(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fix the nan problem above and print real and imaginary part of $M$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sqrt(M).real"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sqrt(M).imag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Matrix multiplication\n",
    "\n",
    "We use the `@` operator to do matrix multiplication with NumPy arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M @ M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute $2I + 3A - AB$ for\n",
    "\n",
    "$$\n",
    "A = \\begin{bmatrix}\n",
    "1 & 3 \\\\\\\n",
    "-1 & 7\n",
    "\\end{bmatrix}\n",
    "\\ \\ \\ \\\n",
    "B = \\begin{bmatrix}\n",
    "5 & 2 \\\\\\\n",
    "1 & 2\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "and $I$ is the identity matrix of size 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1,3],[-1,7]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "B = np.array([[5,2],[1,2]])\n",
    "print(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I = np.eye(2)\n",
    "print(I)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "2*I + 3*A - A@B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Matrix powers\n",
    "\n",
    "There is no symbol for matrix powers and so we must import the function `matrix_power` from the subpackage `numpy.linalg`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.linalg import matrix_power as mpow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = np.array([[3,4],[-1,5]])\n",
    "print(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpow(M,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpow(M,5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare with the matrix multiplcation operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M @ M @ M @ M @ M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpow(M,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M @ M @ M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  2.4 Tranpose\n",
    "\n",
    "We can take the transpose with `.T` attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(M.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that $M M^T$ is a symmetric matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M @ M.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.5 Hermitian conjugate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pauli matrices for spin-$\\frac{1}{2}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigmax = np.array([[0,1],[1,0]])\n",
    "print(sigmax)\n",
    "sigmay = np.array([[0,-1j],[1j,0]])\n",
    "print(sigmay)\n",
    "sigmaz = np.array([[1,0],[0,-1]])\n",
    "print(sigmaz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since Pauli matrices represent physical observable, they must be identical to their Hermitian conjugate $\\hat{\\sigma}_\\alpha=\\hat{\\sigma}_\\alpha^\\dagger$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigmax-np.conjugate(sigmax.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigmay-np.conjugate(sigmay.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigmaz-np.conjugate(sigmaz.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.5 Inverse\n",
    "\n",
    "We can find the inverse using the function `scipy.linalg.inv`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1,2],[3,4]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "la.inv(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.6 Trace\n",
    "\n",
    "We can find the trace of a matrix using the function `numpy.trace`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.trace(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.6 Norm\n",
    "\n",
    "[`scipy.linalg.norm`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.norm.html) returns seven different matrix norms, or one of an infinite number of vector norms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "la.norm(A) # 2-norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "la.norm(A,'fro') # Frobenius norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.7 Determinant\n",
    "\n",
    "We find the determinant using the function `scipy.linalg.det`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1,2],[3,4]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "la.det(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Examples from Mathematics\n",
    "\n",
    "### 3.1 Characteristic polynomials and Cayley-Hamilton theorem\n",
    "\n",
    "The characteristic polynomial of a 2 by 2 square matrix $A$ is\n",
    "\n",
    "$$\n",
    "p_A(\\lambda) = \\det(A - \\lambda I) = \\lambda^2 - \\mathrm{tr}(A) \\lambda + \\mathrm{det}(A)\n",
    "$$\n",
    "\n",
    "The [Cayley-Hamilton Theorem](https://en.wikipedia.org/wiki/Cayley%E2%80%93Hamilton_theorem) states that any square matrix satisfies its characteristic polynomial. For a matrix $A$ of size 2, this means that\n",
    "\n",
    "$$\n",
    "p_A(A) = A^2 - \\mathrm{tr}(A) A + \\mathrm{det}(A) I = 0\n",
    "$$\n",
    "\n",
    "Let's verify the Cayley-Hamilton Theorem for a few different matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trace_A = np.trace(A)\n",
    "det_A = la.det(A)\n",
    "I = np.eye(2)\n",
    "A @ A - trace_A * A + det_A * I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us do this again for some random matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = np.random.randint(0,10,[2,2])\n",
    "print(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trace_N = np.trace(N)\n",
    "det_N = la.det(N)\n",
    "I = np.eye(2)\n",
    "N @ N - trace_N * N + det_N * I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Projection of a vector\n",
    "\n",
    "The formula to project a vector $v$ onto a vector $w$ is\n",
    "\n",
    "$$\n",
    "\\mathrm{proj}_w(v) = \\frac{v \\cdot w}{w \\cdot w} w\n",
    "$$\n",
    "\n",
    "Let's write a function called `proj` which computes the projection $v$ onto $w$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def proj(v,w):\n",
    "    '''Project vector v onto w.'''\n",
    "    v = np.array(v)\n",
    "    w = np.array(w)\n",
    "    return np.sum(v * w)/np.sum(w * w) * w   # or (v @ w)/(w @ w) * w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proj([1,2,3],[1,1,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Linear Systems of equations and alternatives to `scipy.linalg.inv`\n",
    "\n",
    "A [linear system of equations](https://en.wikipedia.org/wiki/System_of_linear_equations) is a collection of linear equations\n",
    "\n",
    "\\begin{align}\n",
    "a_{0,0}x_0 + a_{0,1}x_2 + \\cdots + a_{0,n}x_n & = b_0 \\\\\\\n",
    "a_{1,0}x_0 + a_{1,1}x_2 + \\cdots + a_{1,n}x_n & = b_1 \\\\\\\n",
    "& \\vdots \\\\\\\n",
    "a_{m,0}x_0 + a_{m,1}x_2 + \\cdots + a_{m,n}x_n & = b_m \\\\\\\n",
    "\\end{align}\n",
    "\n",
    "In matrix notation, a linear system is $A \\mathbf{x}= \\mathbf{b}$ where\n",
    "\n",
    "$$\n",
    "A = \\begin{bmatrix}\n",
    "a_{0,0} & a_{0,1} & \\cdots & a_{0,n} \\\\\\\n",
    "a_{1,0} & a_{1,1} & \\cdots & a_{1,n} \\\\\\\n",
    "\\vdots & & & \\vdots \\\\\\\n",
    "a_{m,0} & a_{m,1} & \\cdots & a_{m,n} \\\\\\\n",
    "\\end{bmatrix}\n",
    " \\ \\ , \\ \\\n",
    "\\mathbf{x} = \\begin{bmatrix}\n",
    "x_0 \\\\\\ x_1 \\\\\\ \\vdots \\\\\\ x_n\n",
    "\\end{bmatrix}\n",
    " \\ \\ , \\ \\\n",
    "\\mathbf{b} = \\begin{bmatrix}\n",
    "b_0 \\\\\\ b_1 \\\\\\ \\vdots \\\\\\ b_m\n",
    "\\end{bmatrix} \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 `scipy.linalg.solve`\n",
    "\n",
    "We are mostly interested in linear systems $A \\mathbf{x} = \\mathbf{b}$ where there is a unique solution $\\mathbf{x}$. This is the case when $A$ is a square matrix ($m=n$) and $\\mathrm{det}(A) \\not= 0$. To solve such a system, we can use the function [`scipy.linalg.solve`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html).\n",
    "\n",
    "The function returns a solution of the system of equations $A \\mathbf{x} = \\mathbf{b}$. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1,1],[1,-1]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b1 = np.array([2,0])\n",
    "print(b1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And solve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1 = la.solve(A,b1)\n",
    "print(x1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the output $\\mathbf{x}$ is returned as a 1D NumPy array when the vector $\\mathbf{b}$ (the right hand side) is entered as a 1D NumPy array. If we input $\\mathbf{b}$ as a 2D NumPy array, then the output is a 2D NumPy array. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1,1],[1,-1]])\n",
    "b3 = np.array([[2,2],[0,1]])\n",
    "x3 = la.solve(A,b3)\n",
    "print(x3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 Simple example of two equation with two unknown variables\n",
    "\n",
    "Let's compute the solution of the system of equations\n",
    "\n",
    "\\begin{align}\n",
    "2x + y &= 1 \\\\\\\n",
    "x + y &= 1\n",
    "\\end{align}\n",
    "\n",
    "Create the matrix of coefficients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[2,1],[1,1]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the vector $\\mathbf{b}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.array([1,-1]).reshape(2,1)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And solve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = la.solve(A,b)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can verify the solution by computing the inverse of $A$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ainv = la.inv(A)\n",
    "print(Ainv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And multiply $A^{-1} \\mathbf{b}$ to solve for $\\mathbf{x}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = Ainv @ b\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get the same result. Success!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3 Inverse or solve?\n",
    "\n",
    "It is a bad idea to use the inverse $A^{-1}$ to solve $A \\mathbf{x} = \\mathbf{b}$ if $A$ is large. It is too computationally expensive. Let us create a large random matrix $A$ and vector $\\mathbf{b}$ and compute the solution $\\mathbf{x}$ in 2 ways:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "A = np.random.rand(N,N)\n",
    "b = np.random.rand(N,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the first entries $A$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A[:3,:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And for $\\mathbf{b}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b[:4,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we compare the speed of `scipy.linalg.solve` with `scipy.linalg.inv`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "x = la.solve(A,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "x = la.inv(A) @ b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving with `scipy.linalg.solve` is about twice as fast!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Eigenvalues and eigenvectors\n",
    "\n",
    "Let $A$ be a square matrix. A non-zero vector $\\mathbf{v}$ is an [eigenvector](https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors) for $A$ with [eigenvalue](https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors) $\\lambda$ if\n",
    "\n",
    "$$\n",
    "A\\mathbf{v} = \\lambda \\mathbf{v}\n",
    "$$\n",
    "\n",
    "Rearranging the equation, we see that $\\mathbf{v}$ is a solution of the homogeneous system of equations\n",
    "\n",
    "$$\n",
    "\\left( A - \\lambda I \\right) \\mathbf{v} = \\mathbf{0}\n",
    "$$\n",
    "\n",
    "where $I$ is the identity matrix of size $n$. Non-trivial solutions exist only if the matrix $A - \\lambda I$ is singular which means $\\mathrm{det}(A - \\lambda I) = 0$. Therefore eigenvalues of $A$ are roots of the [characteristic polynomial](https://en.wikipedia.org/wiki/Characteristic_polynomial)\n",
    "\n",
    "$$\n",
    "p(\\lambda) = \\mathrm{det}(A - \\lambda I)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.1 scipy.linalg.eig\n",
    "\n",
    "The function [`scipy.linalg.eig`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eig.html) computes eigenvalues and eigenvectors of a square matrix $A$.\n",
    "\n",
    "Let's consider a simple example with a diagonal matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1,0],[0,-2]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `la.eig` returns a tuple `(eigvals,eigvecs)` where `eigvals` is a 1D NumPy array of complex numbers giving the eigenvalues of $A$, and `eigvecs` is a 2D NumPy array with the corresponding eigenvectors in the columns:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = la.eig(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The eigenvalues of $A$ are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The corresponding eigenvectors are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can [unpack the tuple](https://www.geeksforgeeks.org/unpacking-a-tuple-in-python/):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigvals, eigvecs = la.eig(A)\n",
    "print(eigvals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(eigvecs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we know that the eigenvalues are real numbers (ie. if $A$ is symmetric), then we can use the NumPy array method `.real` to convert the array of eigenvalues to real numbers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigvals = eigvals.real\n",
    "print(eigvals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the position of an eigenvalue in the array `eigvals` correspond to the column in `eigvecs` with its eigenvector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lambda1 = eigvals[1]\n",
    "print(lambda1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v1 = eigvecs[:,1].reshape(2,1)\n",
    "print(v1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A @ v1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lambda1 * v1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2 Symmetric matrices\n",
    "\n",
    "The eigenvalues of a [symmetric matrix](https://en.wikipedia.org/wiki/Symmetric_matrix) are always real and the eigenvectors are always orthogonal! Let's verify these facts with some random matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 4\n",
    "P = np.random.randint(0,10,(n,n))\n",
    "print(P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the symmetric matrix $S = (P + P^T)/2$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = (P + P.T)/2\n",
    "print(S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's unpack the eigenvalues and eigenvectors of $S$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evals, evecs = la.eig(S)\n",
    "print(evals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The eigenvalues all have zero imaginary part and so they are indeed real numbers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evals = evals.real\n",
    "print(evals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The corresponding eigenvectors of $S$ are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(evecs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us check that the eigenvectors are orthogonal to each other:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v1 = evecs[:,0] # First column is the first eigenvector\n",
    "print(v1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v2 = evecs[:,1] # Second column is the second eigenvector\n",
    "print(v2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dot product of eigenvectors $\\mathbf{v}_1$ and $\\mathbf{v}_2$ is zero (the number below is *very* close to zero and is due to rounding errors in the computations), so they are orthogonal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v1 @ v2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.3 Matrix diagonalization\n",
    "\n",
    "A square matrix $M$ is [diagonalizable](https://www.springer.com/gp/book/9783319011943) if it is similar to a diagonal matrix. In other words, $M$ is diagonalizable if there exists an invertible matrix $P$ such that $D = P^{-1}MP$ is a diagonal matrix.\n",
    "\n",
    "A beautiful result in linear algebra is that a square matrix $M$ of size $n$ is diagonalizable if and only if $M$ has $n$ independent eigevectors. Furthermore, $M = PDP^{-1}$ where the columns of $P$ are the eigenvectors of $M$ and $D$ has corresponding eigenvalues along the diagonal.\n",
    "\n",
    "Let's use this to construct a matrix with given eigenvalues $\\lambda_1 = 3, \\lambda_2 = 1$, and eigenvectors $v_1 = [1,1]^T, v_2 = [1,-1]^T$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = np.array([[1,1],[1,-1]])\n",
    "print(P)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = np.diag((3,1)) # Put eigenvectors on diagonal of a matrix\n",
    "print(D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = P @ D @ la.inv(P)\n",
    "print(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evals, evecs = la.eig(M)\n",
    "print(evals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(evecs) # why are eigenvectors not identical to columns of matrix $P$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.4 Matrix powers\n",
    "\n",
    "Let $M$ be a square matrix. Computing powers of $M$ by matrix multiplication\n",
    "\n",
    "$$\n",
    "M^k = \\underbrace{M M \\cdots M}_k\n",
    "$$\n",
    "\n",
    "is computationally expensive. Instead, let's use diagonalization to compute $M^k$ more efficiently\n",
    "\n",
    "$$\n",
    "M^k = \\left( P D P^{-1} \\right)^k  = \\underbrace{P D P^{-1} P D P^{-1} \\cdots P D P^{-1}}_k = P D^k P^{-1}\n",
    "$$\n",
    "\n",
    "Let's compute $M^{20}$ both ways and compare execution time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Pinv = la.inv(P)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "result = M.copy()\n",
    "for _ in range(1,k):\n",
    "    result = result @ M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "P @ D**k @ Pinv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "mpow(M,20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Diagonalization computes $M^{k}$ much faster and, furthermore, it is faster than even built in [`matrix_power`](https://het.as.utexas.edu/HET/Software/Numpy/reference/generated/numpy.linalg.matrix_power.html#numpy.linalg.matrix_power) function!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Examples from Quantum Mechanics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.1 Inner (or dot) product of vectors\n",
    "\n",
    "For dot product of vectors in classical physics, whose components are real numbers, we can use `np.dot` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1,2,3])\n",
    "b = np.array([4,5,6])\n",
    "dot = np.dot(a,b)\n",
    "print(dot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find angle between vectors $\\vec{a}$ and $\\vec{b}$ using dot product and [`scipy.linalg.norm`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.norm.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a_modulus = la.norm(a)\n",
    "b_modulus = la.norm(b)\n",
    "angle = np.arccos(dot/(a_modulus*b_modulus))\n",
    "print(angle*180/np.pi) # angle in degrees"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In quantum mechanics we need dot product of bra (column vector) and ket (row vector), $\\langle v | w \\rangle$, whose components in some representation are complex numbers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wket = np.array([[1,2j,3]]).T\n",
    "wket = wket/la.norm(wket) # normalize\n",
    "print(wket)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vket = np.array([[4,5j,6]]).T\n",
    "vket = vket/la.norm(vket) # normalize\n",
    "vbra = np.conjugate(vket.T)\n",
    "print(vbra)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "overlap = vbra @ wket"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "overlap = overlap.item(0) # convert an array of size 1 to its scalar equivalent\n",
    "print(overlap)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.2 Outer product of vectors yields projection operator <a id='5.2'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Projection operators is given by the outer product of vector with itself, $| a \\rangle \\langle a |$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.T @ a # the result is incorrectly just a scalar instead of expected 3x3 matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "'a' is 1D array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "so it cannot be properly transposed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(a.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a2d=np.array(a).reshape(1,3) # reshape a into row vector as 2D array\n",
    "a2d.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(a2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(a2d.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a2d.T @ a2d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vket @ vbra # projection operator is 3x3 matrix in this case"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.3 Spectral decomposition\n",
    "\n",
    "Any [normal matrix](https://www.springer.com/gp/book/9783319011943), $AA^\\dagger=A^\\dagger A$, is diagonalizable in the sense of having as many eigenvectors as is its dimension or, equivalently, being expressable in terms of the spectral decomposition:\n",
    "\n",
    "$$A=\\sum_n a_n |a_n \\rangle \\langle a_n|$$\n",
    "\n",
    "Let us construct [spectral decomposition](https://www.springer.com/gp/book/9783319011943) for the Pauli $\\hat{\\sigma}_y$ matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D, V = la.eig(sigmay) # find eigenvalues and eigenvectors\n",
    "e1 = V[:,0].reshape(2,1) # first eigenvector reshaped to be 2D array\n",
    "e2 = V[:,1].reshape(2,1) # second eigenvector reshaped to be 2D array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectral = D[0]*(e1 @ np.conjugate(e1.T))  + D[1]*(e2 @ np.conjugate(e2.T))\n",
    "print(spectral)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "la.norm(sigmay-spectral)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.4 Function of diagonalizable matrix <a id='5.4'></a>\n",
    "\n",
    "A function of diagonalizable matrix is obtained from its spectral decomposition\n",
    "\n",
    "$$ F(A)=\\sum_n F(a_n) |a_n \\rangle \\langle a_n|$$\n",
    "\n",
    "Let us find exponent of the Pauli matrix $i{\\sigma}_y \\phi/2$ multiplied by the imaginary unit and angle $\\phi$\n",
    "\n",
    "$$ \\hat{U}=e^{i \\hat{\\sigma}_y \\phi} $$\n",
    "\n",
    "which is the unitary operator of spin rotations around the *y*-axis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%%timeit\n",
    "U = np.exp(1j*D[0])*(e1 @ np.conjugate(e1.T)) + np.exp(1j*D[1])*(e2 @ np.conjugate(e2.T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Udiag = np.diag(D) # put eigenvalues on the diagoanl of matrix Udiag\n",
    "print(Udiag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computational complexity is lower for equivalent unitary transformation of diagonalized matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%%timeit\n",
    "Ualt = V @ np.exp(1j*Udiag) @ np.conjugate(V.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But this is incorrect:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(Ualt)\n",
    "la.norm(U-Ualt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.5 Kronecker product of matrices\n",
    "\n",
    "Kronecker product of matrices A and B is defined by:\n",
    "\n",
    "$$\n",
    "A = \\begin{bmatrix}\n",
    "a_{0,1}B & a_{0,2}B & \\cdots & a_{0,n}B \\\\\\\n",
    "a_{1,0}B & a_{1,1}B & \\cdots & a_{1,n}B \\\\\\\n",
    "\\vdots & & & \\vdots \\\\\\\n",
    "a_{m,0}B & a_{m,1}B & \\cdots & a_{m,n}B \\\\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Kronecker product of matrices is used in quantum mechanics when combining quantum subsystems into composite system. Exchange interaction of two spin-$\\frac{1}{2}$ is described by the Heisenberg Hamiltonian acting in the Hilbert space $\\mathcal{H}_1^\\mathrm{spin} \\otimes \\mathcal{H}_2^\\mathrm{spin}$\n",
    "\n",
    "$$ \\hat{H} = \\hat{\\vec{\\sigma}}_1 \\cdot \\hat{\\vec{\\sigma}}_2 $$\n",
    "\n",
    "Let us find its matrix representation using the vector of the Pauli matrices, $\\hat{\\vec{\\sigma}}_1 = \\hat{\\vec{\\sigma}} \\otimes \\hat{I}_{2 \\times 2}$ and  $\\hat{\\vec{\\sigma}}_2 = \\hat{I}_{2 \\times 2} \\otimes \\hat{\\vec{\\sigma}}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I2 = np.eye(2,2)\n",
    "np.kron(sigmax,I2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.kron(I2,sigmax) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = np.kron(sigmax,I2) @ np.kron(I2,sigmax) + np.kron(sigmay,I2) @ np.kron(I2,sigmay) + np.kron(sigmaz,I2) @ np.kron(I2,sigmaz)\n",
    "print(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.6 Tight-binding Hamiltonian of a disordered chain of four atoms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use Anderson disorder as the uniform random variable $\\varepsilon_i \\in [-W/2,W/2]$ on the diagonal of tight-binding Hamiltonian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = 4.\n",
    "H=np.diag(W*(np.random.rand(4)-1/2),k=0)\n",
    "print(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add hoppings between the nearest-neighbor sites:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = H + np.diag(-np.ones(3), k=-1)  # add hopping on the lower diagonal\n",
    "H = H + np.diag(-np.ones(3), k=1)   # add hopping on the upper diagonal\n",
    "print(H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "la.norm(H-np.conjugate(H.T)) # Hamiltonian is always Hermitian matrix, so this check can prevent bugs in the code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Send corrections to <code>bnikolic@udel.edu</code>"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}