{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 6\n",
    "\n",
    "* Built-in Linear Algebra routines.\n",
    "* Vector-vector, matrix-vector products."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q & A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# difference between the following two functions\n",
    "def myfunc1(x):\n",
    "    out = 0\n",
    "    for i in range(len(x)):\n",
    "        out += x\n",
    "        return out\n",
    "    \n",
    "def myfunc2(x):\n",
    "    out = 0\n",
    "    for i in range(len(x)):\n",
    "        out += x\n",
    "    return out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vector-vector inner product\n",
    "The inner product (or dot product): for $\\mathbf{a}, \\mathbf{b} \\in \\mathbb{R}^n$, their inner product is\n",
    "$\\mathbf{a}\\cdot\\mathbf{b} = \\displaystyle\\sum_{i=1}^n a_i b_i.$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for loop version\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# semi-vectorized numpy version\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# built-in function\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use function built-in the ndarray object\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matrix vector multiplication\n",
    "We have $\\mathbf{A}\\in \\mathbb{R}^{m\\times n}$ being a matrix, and $\\mathbf{x}\\in \\mathbb{R}^n$ being a (column) vector. Then we can compute $\\mathbf{A}\\mathbf{x}$:\n",
    "$$\\mathbf {A} \\mathbf{x}=\n",
    "\\begin{pmatrix}a_{11}&a_{12}&\\cdots &a_{1n}\\\\a_{21}&a_{22}&\\cdots &a_{2n}\n",
    "\\\\\\vdots &\\vdots &\\ddots &\\vdots \\\\a_{m1}&a_{m2}&\\cdots &a_{mn}\\\\\\end{pmatrix}\n",
    "\\begin{pmatrix}x_{1}\\\\x_{2}\\\\\\vdots \\\\x_{n}\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}a_{11}x_{1}+\\cdots +a_{1n}x_{n}\\\\a_{21}x_{1}+\\cdots +a_{2n}x_{n}\\\\ \\vdots \\\\a_{m1}x_{1}+\\cdots +a_{mn}x_{n}\\end{pmatrix}\n",
    "= \\begin{pmatrix} \\mathbf{a}_1\\cdot \\mathbf{x} \\\\ \\mathbf{a}_2\\cdot \\mathbf{x} \\\\ \\vdots \\\\ \n",
    "\\mathbf{a}_m\\cdot \\mathbf{x}\\end{pmatrix}\n",
    "$$\n",
    "where the matrix $\\mathbf {A} $ has the following representation:\n",
    "$$\n",
    "\\mathbf{A} = (\\mathbf{a}_1 , \\mathbf{a}_2 ,\\cdots ,\\mathbf{a}_m)^{\\top}.\n",
    "$$\n",
    "with $\\mathbf{a}_j$ being the column vector representing $\\mathbf{A}$'s $j$-th row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1,2,3],[4,5,6],[7,8,9]])\n",
    "x = np.array([1, 0, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ???\n",
    "A*x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.multiply(A,x) # same with A*x, multiply arrays element-wise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# more example\n",
    "S = np.arange(9).reshape(3,3)\n",
    "v = np.arange(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.multiply(S,v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Question: how to we implement the matrix vector multiplication?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for loop version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# built-in function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1: Convolution in 2D using a filter\n",
    "\n",
    "Suppose we have a big matrix `A` and a small filter matrix (called Kernel) `K`:\n",
    "<img src=\"convolution1.JPG\" alt=\"Convolution\" width=\"500\"/>\n",
    "It is adding up the element-wise product of every 3 by 3 block in the big matrix `A` with the kernel `K`:\n",
    "<img src=\"convolution2.JPG\" alt=\"Convolution\" width=\"700\"/>\n",
    "\n",
    "* Implement the a convolution using a `3x3` filter `K` with a `(28,20)` array `A` that are generated using the following code:\n",
    "```python\n",
    "A = np.random.randint(90,110, size=(28,20))\n",
    "K = np.array([[0,-1,0], [-1,5,-1], [0,-1,0]])\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving linear system\n",
    "Now let us consider when $m=n$, we want to solve a system of linear equations:\n",
    "$$\\mathbf{A} \\mathbf{x} = \\mathbf{b},$$\n",
    "which is\n",
    "$$\\left\\{\n",
    "\\begin{aligned} a_{11}x_{1}+a_{12}x_{2}+\\cdots +a_{1n}x_{n}&=b_{1}\\\\a_{21}x_{1}+a_{22}x_{2}+\\cdots +a_{2n}x_{n}&=b_{2}\\\\&\\ \\ \\vdots \\\\a_{n1}x_{1}+a_{n2}x_{2}+\\cdots +a_{nn}x_{n}&=b_{n}.\\end{aligned}\n",
    "\\right.\n",
    "$$\n",
    "Instead of manual implementing the pivoting process using `for` loops, we can use the Linear Algebra routines `numpy.linalg` in NumPy, and later using `scipy.linalg` in SciPy (a comprehensive scientific computing library).\n",
    "\n",
    "<br><br>\n",
    "Reference: [https://docs.scipy.org/doc/numpy-1.15.1/reference/routines.linalg.html](https://docs.scipy.org/doc/numpy-1.15.1/reference/routines.linalg.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import linalg as LA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first we have to check the rank\n",
    "A = np.array([[1,2,3],[4,5,6],[7,8,9]])\n",
    "LA.matrix_rank(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[2,-1,1],[1,2,-1],[0,-1,2]])\n",
    "LA.matrix_rank(A) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.array([1,1,1])\n",
    "x = LA.solve(A,b)\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other common function in linear algebra\n",
    "The determinant, inverse, and eigenvalues/eigenvectors: `det`, `inv`, `eig`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ainv = LA.inv(A)\n",
    "print(Ainv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LA.det(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w, v = LA.eig(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2: Eigenvectors of tri-diagonal matrices\n",
    "\n",
    "The matrix $A$ is as follows:\n",
    "$$A = \n",
    "\\begin{pmatrix} 2 & -1 \n",
    "\\\\ -1 & 2 &-1 \n",
    "\\\\ & -1 & 2 & -1\n",
    "\\\\ && \\ddots & \\ddots & \\ddots\n",
    "\\\\ &&& -1 & 2 & -1\n",
    "\\\\& &&&  -1 & 2\n",
    "\\end{pmatrix}$$\n",
    "* Use `numpy` function `eye` to generate this matrix without any loops (reference: [numpy.eye](https://docs.scipy.org/doc/numpy/reference/generated/numpy.eye.html)).\n",
    "* Use the `eig` function in `linalg` to find the eigenvalues `eigvals` and the eigenvectors `eigvecs` of this matrix.\n",
    "* Generate 32 equi-distant points on $[0,\\pi]$ using `linspace`, plot the eigenvectors which correspond to the smallest three eigenvalues against these points using `plot` function in `matplotlib.pyplot` (reference: Lecture 5), what does these eigenvectors look like?"
   ]
  }
 ],
 "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.7.2"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}