14_Basics_of_Linear_Algebra.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "14_Basics of Linear Algebra.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/fmml2021/fmml2021/blob/main/14_Basics_of_Linear_Algebra.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZfkNibSQgtCZ"
      },
      "source": [
        "## Basics of Linear Algebra"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bAPP654PeaPa"
      },
      "source": [
        "Import the NumPy library for Linear Algebra functions and Matplotlib for some plotting functions. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZHOzQlD7eQxm"
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YXhIDS8NevLp"
      },
      "source": [
        "### Transpose of Matrix\n",
        "\n",
        "Matrix transpose is performed with the transpose method on a nested list or a Python array, or a higher-dimensional Numpy array. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NbJbiJ6Sewe0"
      },
      "source": [
        " # Transpose of a Matrix (as nested list)\n",
        " a = [[1,2,3,4],[2,3,4,5]]\n",
        " b = np.transpose(a)\n",
        " print('a\\n',a)\n",
        " print('b\\n',b) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3fUxqcc7fDHn"
      },
      "source": [
        "If the matrix is a NumPy array, it can be treated as an object and method T can be applied over it as follows."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9lWWJwbsfEKj"
      },
      "source": [
        "  # Transpose of a Matrix (as NumPy array)\n",
        "print (\"Matrix and its Transpose\")\n",
        "a = np.array([[1,2,3,4],[2,3,4,5]])\n",
        "b = a.T\n",
        "print('a\\n',a)\n",
        "print('b\\n',b) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YFEGttjsfiQe"
      },
      "source": [
        "The dot method of NumPy performs dot-matrix product (scalar product) for 1D or higher dimensional arrays. If the inputs are scalars (numbers), it performs multiplication."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "euFPFMmNfjG-"
      },
      "source": [
        " # scalars\n",
        " a = 5\n",
        " b = 3\n",
        " z = np.dot(a,b)\n",
        " print(z) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bXk9XEnqfuF4"
      },
      "source": [
        " z = a * b\n",
        " print(z) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kqJ5cmN6f4Go"
      },
      "source": [
        "In the case of one- or higher-dimensional arrays, the inputs can be either NumPy arrays, Python arrays, Python lists or Python’s nested lists."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "U2dUeiHLf5Bx"
      },
      "source": [
        " # 1D arrays or vectors\n",
        " a = np.array([1,2,3])  # or a = [1,2,3]\n",
        " b = np.array([2,3,4])  # or b = [2,3,4]\n",
        " z = np.dot(a,b)\n",
        " print(z) \n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kWFmevo-gIWV"
      },
      "source": [
        " # 2D arrays or matrices\n",
        " a = [[1,2,3],[2,0,3],[7,-5,1]]\n",
        " b = [[3,-1,5],[-2,-6,4], [0,4,4]]\n",
        " z = np.dot(a,b)\n",
        " print(z) \n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_cEG72RdIkmI"
      },
      "source": [
        "We can obtained the same result using `np.matmul()`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hN3aQ5jVIais"
      },
      "source": [
        "z = np.matmul(a,b)\n",
        "print(z)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QOgLFDj8g31t"
      },
      "source": [
        "### Numpy Arrays\n",
        "A NumPy array is a Numpy object upon which the dot method can be performed as below. However, this method accepts only NumPy arrays to operate on.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "D61xhcT_g7Eu"
      },
      "source": [
        " # convert lists into NumPy arrays\n",
        " newa = np.array(a)\n",
        " newb = np.array(b)\n",
        " z = newa.dot(newb)\n",
        " print(z) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QQAta3Razf2b"
      },
      "source": [
        "### The multi_dot method\n",
        "\n",
        "It performs dot (scalar) product with 2 or more input matrices. First and last arrays can be either 1D or 2D arrays. However, the dimensions of the matrices must suit subsequent scalar matrix multiplication. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PON5zE_2zXaT"
      },
      "source": [
        " # matrices with random integers: entries ranging from -4 to 4\n",
        " a = np.random.randint(-4,4,(500,5))\n",
        " b = np.random.randint(-4,4,(5,1000))\n",
        " c = np.random.randint(-4,4,(1000,10))\n",
        " d = np.random.randint(-4,4,(10,2000))\n",
        " e = np.random.randint(-4,4,(2000,200))\n",
        " # Perform multiple matrix multiplication\n",
        " z = np.linalg.multi_dot([a,b,c,d,e]) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Nb9iak7Pz4cg"
      },
      "source": [
        "The result of this method can be obtained with successive dot products of matrices but multi_dot functions in an optimized manner. It decides the order of dot multiplication to complete the entire process efficiently."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fHMnlw7BzzvO"
      },
      "source": [
        " %%time\n",
        " z = np.linalg.multi_dot([a,b,c,d,e])\n",
        " print(z, '\\n') "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "iBUZFx8_0EDq"
      },
      "source": [
        " %%time\n",
        " z = a.dot(b).dot(c).dot(d).dot(e)\n",
        " print(z, '\\n') \n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "51oLesAq0MnM"
      },
      "source": [
        "#### Observe the CPU time that multi_dot consumes as against CPU time that successive dot methods consume to arrive at the same solution. Which one is more efficient, why ?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uCefpm7o0j5y"
      },
      "source": [
        "### Inner Product\n",
        "\n",
        "The inner product is the scalar multiplication of one vector (or matrix) and the transpose of another vector (or matrix). If both arrays are 1D arrays, their dimensions should be identical. If either or both arrays are higher-dimensional, then the last dimensions of both arrays should be identical."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bmLy4CCs0gCl"
      },
      "source": [
        " a = np.array([[1,2,3], [4,-1,0]])\n",
        " b = np.array([6,3,2])\n",
        " z = np.inner(a,b)\n",
        " print(z) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0aVTzlXC1AwG"
      },
      "source": [
        "\n",
        "# The same results can be obtained using the dot method as follows.\n",
        "\n",
        "a.dot(b.T)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KX34ZLjv1JHD"
      },
      "source": [
        "### Outer Product\n",
        "\n",
        "Outer product is the dot product of a column vector of size Mx1 and a row vector of size 1xN. The resulting array is a matrix of size MxN.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hktyza1R1G0K"
      },
      "source": [
        " a = np.array([1,2,3,4,5])\n",
        " b = np.array([6,3,2])\n",
        " z = np.outer(a,b)\n",
        " print(z) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5vsxEYGq1pKl"
      },
      "source": [
        "The last dimension of the second array and the second-to-last dimension of the first array should be identical to perform matrix multiplication. Further, the symbol @ is also used to perform matrix multiplication."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ltAaQjPt1p95"
      },
      "source": [
        " # Here, 'a' matrix is 3D, which means there are 3 matrices each of 2x5 size\n",
        " # Similarly, for 'b' matrix\n",
        " # So we perform 3 matrix multiplication operations each with 2x5 and 5x3 matrices from a and b\n",
        " a = np.random.random([3,2,5])\n",
        " b = np.random.random([3,5,3])\n",
        " z = a @ b\n",
        " print(z.shape) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tKi7vaCc1xXd"
      },
      "source": [
        "### Matrix Determinant\n",
        "\n",
        "Matrix determinant can be calculated using the method det available in the linalg module."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ONGBxJuS1zGP"
      },
      "source": [
        " # generate a random integer matrix of size 3 by 3\n",
        " a = np.random.randint(1,10,[3,3])\n",
        " det = np.linalg.det(a)\n",
        " print(int(det)) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kmX7izKQ18yY"
      },
      "source": [
        "### Matrix Inverse\n",
        "\n",
        "Inverse of a square matrix can be derived using the inv method of the linalg module.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bPvi1QGq1-g8"
      },
      "source": [
        " a = np.random.randint(1,10,[3,3])\n",
        " inv = np.linalg.inv(a)\n",
        " print(a)\n",
        " print()\n",
        " print(inv) \n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2w2sxTzF2Dpw"
      },
      "source": [
        "### Matrix Power\n",
        "\n",
        "Matrix Power is a general method to obtain either positive or negative powers of a given square matrix. The first negative power of a matrix is technically termed its inverse. Thus, the matrix_power method can be used to find the inverse or any power of a matrix."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "uHbQXzDd2JDJ"
      },
      "source": [
        " a = np.random.random([4,4])\n",
        " # positive powers of matrix\n",
        " a_2 = np.linalg.matrix_power(a, 2)\n",
        " a_7 = np.linalg.matrix_power(a, 7)\n",
        " # inverse of matrix\n",
        " a_inv_1 = np.linalg.matrix_power(a, -1)\n",
        " a_inv_3 = np.linalg.matrix_power(a, -3)\n",
        " print('matrix \\n', a)\n",
        " print('\\n matrix to the power 2\\n', a_2)\n",
        " print('\\n matrix to the power 7\\n', a_7)\n",
        " print('\\n matrix inverse \\n', a_inv_1)\n",
        " print('\\n matrix cubic inverse \\n', a_inv_3) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JoJmR88k2ObP"
      },
      "source": [
        "### Eigenvalues and Eigenvectors\n",
        "\n",
        "Eigenvalues and Eigenvectors of a matrix can be determined as follows. If Eigen values cannot be determined, the method throws an error (Eg. Singular matrix)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ye8paxl82P0x"
      },
      "source": [
        " a = np.arange(9).reshape(3,3)\n",
        " eig_val, eig_vec = np.linalg.eig(a)\n",
        " print('Eigenvalues are: \\n', eig_val)\n",
        " print('\\nEigenvectors are: \\n', eig_vec) \n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "iyMHkixI2X-Y"
      },
      "source": [
        "# Eigenvalues alone can be determined using the method eigvals as shown below.\n",
        "\n",
        "a = np.arange(9).reshape(3,3)\n",
        "eigenvalues = np.linalg.eigvals(a)\n",
        "print(eigenvalues) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9ftxHuxC2eSE"
      },
      "source": [
        "### Traces of a Matrix\n",
        "\n",
        "Traces of a square matrix is the summation of its diagonal elements."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "WfqE1uaA2fqh"
      },
      "source": [
        "a = np.eye(5)\n",
        "print(a)\n",
        "z = np.trace(a)\n",
        "print('\\nTrace of matrix is: ',z) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8I1VmBae3QSQ"
      },
      "source": [
        "### Matrix Norm\n",
        "\n",
        "Matrix or vector norm is calculated using the norm method of the linalg module.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "p1FgJ1683RJ-"
      },
      "source": [
        "a = np.arange(12).reshape(4,3)\n",
        "z = np.linalg.norm(a)\n",
        "print(a)\n",
        "print('\\n Frobenius Norm of above matrix:')\n",
        "print(z) \n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Wp8kegE03Y1U"
      },
      "source": [
        "#### Norm of Matrix\n",
        "\n",
        "Axis-wise norm determination is also possible by specifying the axis as an integer.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-NDr7tMa3aKu"
      },
      "source": [
        "# Norm along axis 0\n",
        "a = np.arange(12).reshape(4,3)\n",
        "z = np.linalg.norm(a, axis=0)\n",
        "print(z) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rPOUS32B3q-c"
      },
      "source": [
        "### Solving System of Equations\n",
        "\n",
        "When we think of Linear Algebra, the system of linear equations comes to our mind first, as it is tedious, time-consuming and error-prone. NumPy solves systems of linear equations in a fraction of seconds!"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bM2EczhP3saM"
      },
      "source": [
        "# Coefficient Matrix\n",
        "a = np.random.randint(1,20,[4,4])\n",
        "# Dependent variable vector\n",
        "b = np.array([4,9,12,7])\n",
        "# solution\n",
        "x = np.linalg.solve(a,b)\n",
        "print('Coefficient Matrix')\n",
        "print(a)\n",
        "print('\\nDependent Variable vector')\n",
        "print(b)\n",
        "print('\\nSolution')\n",
        "print(x) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-b-DT76N4B_Q"
      },
      "source": [
        "# Check for correctness\n",
        "B = a.dot(x)\n",
        "print(B) \n",
        "\n",
        "# This ‘B’ array is identical to the input ‘b’ array. Hence, our solution is correct. "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "N9A11eBh4K6F"
      },
      "source": [
        "### Singular Value Decomposition\n",
        "\n",
        "Singular Value Decomposition (SVD) is one of the great dimension-reduction algorithms in machine learning. It identifies the principal components and arranges them by rank. The top ranked components contribute greatly to the original array. Here, we explore SVD with an image to get better understanding."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4bcLK0uw4MSF"
      },
      "source": [
        "from skimage import data\n",
        "# download a sample image\n",
        "image = data.astronaut()\n",
        "print(image.dtype, image.min(), image.max(), image.shape)\n",
        "plt.imshow(image)\n",
        "plt.show() "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jmTyUD7F4Kpg"
      },
      "source": [
        "# Normalize the image by dividing it by the maximum value, 255 and reorder the axes to be (3,400,600).\n",
        "\n",
        "# normalize image\n",
        "img = image/255.0\n",
        "# reorder the axes to have colour channels as the first dimension\n",
        "img = np.transpose(img, axes=(2,0,1)) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OSht9dHH45WB"
      },
      "source": [
        "Perform SVD on the image. It decomposes the original image into three components: U matrix, Sigma vector, and V matrix. The Sigma vector is the diagonal entries of the Sigma matrix. Hence, it is advisable that the Sigma vector may be reformed into a diagonal matrix. It should be noted that the first dimension 3 refers to the three colour channels."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mpHJe2Pk46IA"
      },
      "source": [
        "U,S,V = np.linalg.svd(img)\n",
        "print(U.shape, S.shape, V.shape) \n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "656YnLf94Kn7"
      },
      "source": [
        "# S matrix should have dimensions suitable for matrix multiplication\n",
        "Sigma = np.zeros((3,512,512))\n",
        "for i in range(3):\n",
        "  np.fill_diagonal(Sigma[i,:,:], S[i,:])\n",
        "print(Sigma.shape) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-p87agUP5FjH"
      },
      "source": [
        "# Reconstruct the original image without any dimension reduction.\n",
        "reconst = U @ Sigma @ V\n",
        "reconst = np.transpose(reconst, axes=(1,2,0))\n",
        "plt.imshow(reconst)\n",
        "plt.show() \n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FrINdxMR5eOZ"
      },
      "source": [
        "# NumPy SVD reconstruction\n",
        "\n",
        "# Reconstruct the data by reducing the common dimensions from 512 to 50.\n",
        "\n",
        "k = 50\n",
        "reconst = U @ Sigma[:,:,:k] @ V[:,:k,:]\n",
        "reconst = np.transpose(reconst, axes=(1,2,0))\n",
        "plt.imshow(reconst)\n",
        "plt.show() "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tRD-0vSM6a1u"
      },
      "source": [
        "#### Dimensionality reduction\n",
        "\n",
        "It is amazing that we reconstructed the image with most details, even at a reduction of ¼. \n",
        "\n",
        "We can once again reconstruct the same image by reducing the data points from 512 to 20.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "F3vbxx046iPy"
      },
      "source": [
        "k = 20\n",
        "reconst = U @ Sigma[:,:,:k] @ V[:,:k,:]\n",
        "reconst = np.transpose(reconst, axes=(1,2,0))\n",
        "plt.imshow(reconst)\n",
        "plt.show() "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q1TbCX8p6t9x"
      },
      "source": [
        "#### SVD reconstruction\n",
        "\n",
        "Out of 400 data points, merely 20 data points can reconstruct the image with key feature details! This is why the old mathematical algorithm- SVD is still popular in machine learning."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CHmFf8GfKGK8"
      },
      "source": [
        "### Additional References\n",
        "\n",
        "[Numpy Matrix Multiplication Terminology](https://likegeeks.com/numpy-matrix-multiplication/#Basic_Terminologies)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1A3npCtyK5eT"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}