project#6.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1afd7c0e-02ea-4e93-9c3b-cf54faaffe7f",
   "metadata": {},
   "source": [
    "# Math 266 - Project #6 - Singular value decomposition\n",
    "\n",
    "Execute each cell below and follow the instructions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca30a01f-afc5-4a58-9e15-aa2c5ba13c1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "from numpy import cos,sin,pi\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy.linalg import norm\n",
    "import imageio\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "from sympy import *\n",
    "from sympy.interactive.printing import init_printing\n",
    "init_printing(use_unicode=False, wrap_line=False)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18afe9ab-d9e7-43b0-8ba9-dfa66407288a",
   "metadata": {
    "tags": []
   },
   "source": [
    "---\n",
    "# The singular value decomposition for a matrix A is given by:\n",
    " $$A = U \\Sigma V^T$$\n",
    " \n",
    " Where U and V are unitary matrices. $\\Sigma$ is a \"diagonal matrix containing the singular values for A. The singular values $\\sigma_i$ are in descending order in matrix $\\Sigma$ .\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9449a31-ef1c-4334-b70f-3b7bf16838fb",
   "metadata": {},
   "source": [
    "# Example from text.  Construct a Singular Value Decomposition of the matrix A,  p. 443 in text.\n",
    "\n",
    "$$ A=\\begin{bmatrix}4 & 11 & 14 \\\\8 & 7 & -2 \\end{bmatrix}\\ $$\n",
    "\n",
    "$$A = U \\Sigma V^{T}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e8613b0-c2da-4327-adf7-ee07490253b0",
   "metadata": {
    "tags": []
   },
   "source": [
    "---\n",
    "### To construct the singular value for the matrix A the text outlines the following steps.\n",
    "\n",
    "Step 1. Find the orthogonal diagonalization of $A^{T}A$.\n",
    "\n",
    "step 2. Form matrix V from the eigenvectors for $A^{T}A$.\n",
    "\n",
    "step 3. The singular values are the square root of eigenvalues of $A^{T}A$.<br> &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;\n",
    "&nbsp;&nbsp;&nbsp;&nbsp; Put the values in descending order and form the matrix $\\Sigma$.\n",
    "\n",
    "step 4. Finally form the matrix U.  Column i of U is given by:    $\\vec{\\mathbf{u_i}}=A\\vec{\\mathbf{v_i}}/\\sigma_i$\n",
    "\n",
    "---\n",
    "\n",
    "\n",
    "### Using python we can find the singular value decomposition for matrix A with the simple command:\n",
    "``` python\n",
    "np.linalg.svd(A)\n",
    "\n",
    "```\n",
    "---\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf189dc6-dfd9-45bb-9316-7158d5e4486f",
   "metadata": {},
   "source": [
    "Execute the cell below to create matrix A."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba812465-3edc-4cc1-adb5-b6dd1ac9f7ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([4,11,14,8,7,-2]).reshape(2,3)\n",
    "print (\"A=\",end=\"\");Matrix(A) # command to simply display matrix in nice format"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "797fe821-eb8b-4506-b723-71fc9402c499",
   "metadata": {},
   "source": [
    "___\n",
    "Execute cell below to compute the singular value decomposition of A and view the results. Compare the result to page 444 in text.\n",
    "\n",
    "___"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7cf0734-0240-410e-8027-effb57f4781b",
   "metadata": {},
   "outputs": [],
   "source": [
    "u,s,v = la.svd(A)  # easier way to do it\n",
    "\n",
    "[Matrix(u),Matrix(np.diag(s)),Matrix(v)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4739031-22c7-4446-8bd0-e98527b10b07",
   "metadata": {},
   "source": [
    "Now we will confirm that $A = U \\Sigma V^{T}$ .  Note the svd command returns an array s of singular values and we need to create the diagonal matrix $ \\Sigma $ with the command:\n",
    "```python\n",
    "np.diag(s)\n",
    "\n",
    "```\n",
    "Now $\\Sigma$ must be the same size as A so we need to append a column of zeros to make $\\Sigma$ size (3,2). <br>\n",
    "\n",
    "We do this with the command:\n",
    "\n",
    "```python\n",
    "np.c_[np.diag(s),[0,0]]\n",
    "\n",
    "```\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5724627c-f4ac-43fb-9bdf-af75bc9b6750",
   "metadata": {},
   "source": [
    "Execute the cells below to confirm the decomposition. That is that  $A = U \\Sigma V^{T}$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8a3a990-55be-498c-ab91-27e47a5e8efc",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma=np.c_[np.diag(s),[0,0]]\n",
    "print(\"sigma = \");Matrix(sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ae1f504-073b-4ff3-8146-7bc74673ec16",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(u@sigma@v) # confirm the decomposition - see the output below"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbd8b0bc-0c72-4687-84e5-7857befbd27f",
   "metadata": {},
   "source": [
    "---\n",
    "---\n",
    "\n",
    "# Application 1.    Image compression\n",
    "\n",
    "Recall that a \"grayscale\" image is just an rectangular array (matrix) of grayscale values for each pixel.\n",
    "\n",
    "Th rank n approximation of an Image is given by:\n",
    "\n",
    "$$A = \\sum_{i=0}^{n-1}\\sigma_i\\vec{\\mathbf{u_i}}\\vec{\\mathbf{v_i}}^{T} $$\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce228a65-b292-41b5-8634-9ea8ad34febc",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "Execute the following cells to import the test image and display it. This is an image of the mathematician Bernhare Reimann. We store the image <br>\n",
    "in a (246,225) matrix and call it \"img\". \n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a03287b3-124a-4595-8b10-f13250b587e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "img = imageio.imread(\"https://github.com/rmartin977/Math-266/blob/main/rieman.jpeg?raw=true\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a8b1ddf-4c7a-4b07-98ee-1bffc4ba2de3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "img.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96935a96-413c-48c0-a529-00cc8ac4934d",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(img,cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86a30002-78b1-420e-b4a8-36618510f4f8",
   "metadata": {},
   "source": [
    "---\n",
    "Execute the following to compute SVD of image.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccbf2f7e-3d00-4704-b288-6c9afd399957",
   "metadata": {},
   "outputs": [],
   "source": [
    "U,S,V = la.svd(img)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b0ff3da-7090-44d7-b69b-6dffb3a9a4ce",
   "metadata": {},
   "source": [
    "---\n",
    "Below we do a plot of the magnitudes of the singular values. Notice how they fall off rapidly. This tells us a low rank approximation (image compression) will be good.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbe18976-e67d-446d-941d-374908a4fafb",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(S)\n",
    "plt.grid()\n",
    "plt.title(\"singular values\");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57dde2a0-235f-4055-83fe-612d8aebd36f",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "Execute the cells below to see a rank 1 approximation for our image.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac33221f-f989-43ec-841a-1cdeca839e1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "rank_1 = S[0]*np.outer(U[:,0],V[0,:])\n",
    "rank_1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6bbf3a1-1524-4311-9f27-ddcdea65422d",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(rank_1,cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "497fa5cf-a7df-409f-a4e6-01131ef8244b",
   "metadata": {},
   "source": [
    "Below we define a helper function that will compute a rank_k approximation for our images. The function basically computes the partial sum:\n",
    "\n",
    "$$ \\sum_{i=0}^{k-1}\\sigma_i\\vec{\\mathbf{u_i}}\\vec{\\mathbf{v_i}}^{T} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa8f343a-b781-493b-9d64-d11ae441cbdd",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rank(k):\n",
    "    '''\n",
    "    The function will return a rank_k approximation for the image of Rieman.\n",
    "    '''\n",
    "    return np.sum([ S[i]*np.outer(U[:,i],V[i,:]) for i in range(k)],axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d0d5787-9e16-47a5-907a-e2539cb1f93c",
   "metadata": {},
   "source": [
    "The cell below computes a rank 5 approximation of the image matrix and the result is displayed. Already starting to look like Bernhard. Execute the cell to output.  Change the rank and oberve difference in output.  Try rank = 50."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c294592-1fae-4951-8fd4-f6cccea824ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(rank(5),cmap='gray')  # rank 5  change the value to see different rank approximations."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7717f3a7-a3d3-415e-b8b6-4b6459fa6f83",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "Below we create a rank-5 approximation and store it in matrix Berhnard_5 and save image to file.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5141db12-94d8-4925-99a5-abe26c67a8ea",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "Bernhard_5 = rank(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70f6b781-ffa5-40bb-b76e-7887ad85804d",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.axis('off')\n",
    "plt.imshow(rank_5,cmap='gray')\n",
    "plt.savefig(\"Bernhard_5.jpg\") # you will need to change path to file\n",
    "\n",
    "# to save on your G-drive you need to:\n",
    "# 1. Mount your drive.\n",
    "# 2. Change the path to the file as \"drive/MyDrive/Bernhard_5.jpg\"\n",
    "# 3. Contact me if you need help. Verify that your image is saved and you can display it outside of google colab."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00622843-c6bf-4e8e-b852-465f56299686",
   "metadata": {},
   "source": [
    "---\n",
    "---\n",
    "\n",
    "# Your turn. \n",
    "\n",
    "Import the image of MonaLisa and creat a rank_5 approximation for this image.  Save a copy of the image to your G-drive.  You will need to upload the image to Gradescope later.\n",
    "\n",
    "---\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cf82aae-718b-455e-a4d9-b6c116874eba",
   "metadata": {},
   "outputs": [],
   "source": [
    "image = imageio.imread(\"https://github.com/rmartin977/Math-266/blob/main/ml.jpg?raw=true\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05908510-234d-49e0-a5fc-68f08b7ee397",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(image,cmap='gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f9dce56-e16e-4a75-a621-57aff56e0917",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Enter your code here to create a rank-5 approximation of monalisa and save the jpeg to your G Drive.  Name the file monalisa_5.jpg ."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e2becd3-68a1-4fc1-8f17-cd456809b6f2",
   "metadata": {},
   "source": [
    "---\n",
    "---\n",
    "# Application #2. Pseudo inverse.\n",
    "\n",
    "We will find a least square solution for the system from exam #3. $$ \\begin{bmatrix}1 & 1  \\\\1 & 2\\\\1 & 3\\end{bmatrix}\\ \\begin{bmatrix}\\beta_0  \\\\ \\beta_1  \\end{bmatrix}= \\begin{bmatrix}0  \\\\ 1\\\\1  \\end{bmatrix} $$ \n",
    "\n",
    "By computing the pseudo inverse of matrix A.  A non square matrix does not have a true inverse however we can compute a pseudo inverse by taking the inverse of the \"reduced\" singular value decomposition of A.  This inverse is denoted $A^t$ and is computed as follows.\n",
    "\n",
    "$$A^{t} = V\\Sigma ^{-1} U^{T}$$\n",
    "\n",
    "\n",
    "\n",
    "---\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c63177f3-4b42-4b88-a585-b6114bd20773",
   "metadata": {},
   "source": [
    "---\n",
    "First we recompute the svd for A, and then we comptute the pseudo inverse. We will call the pseduo inverse A_pinv. Execute the cells below.\n",
    "\n",
    "___\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebcdd1ba-6502-4368-b8e4-c975d897438c",
   "metadata": {},
   "source": [
    "Define A and b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "247f25aa-5158-4a81-a3ca-194ea0699a4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([1,1,1,2,1,3]).reshape(3,2)\n",
    "b = np.array([0,1,1])\n",
    "Matrix(A),Matrix(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d06418e-f868-4a95-92bf-1c78ed38fcde",
   "metadata": {},
   "source": [
    "Compute pseudo inverse.  See 7.4 in the text."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df77ab53-2354-4586-b656-27271d4289e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "u,s,v = la.svd(A,full_matrices=False) # compute reduced svd see 7.4 in text\n",
    "sigma_inv = np.diag(1/s)\n",
    "A_pinv = v.T@sigma_inv@u.T\n",
    "Matrix(A_pinv)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec88ef95-a6fd-4d5a-9964-bc84a25b3434",
   "metadata": {},
   "source": [
    "Numpy has a command to compute the pseudo inverse directly.  The command is la.pinv(A).  Execute the cell below and compare."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37b84f80-e2a4-4bc7-8a82-b6173320958b",
   "metadata": {},
   "outputs": [],
   "source": [
    "la.pinv(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "936d886c-b7ee-4700-842e-22f74862a4cc",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "Finally let us use A_inv to solve system. Execute the cell below.\n",
    "\n",
    "___"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc107be1-02bc-4be2-a7a9-bee674dece69",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_hat = A_pinv @ b\n",
    "Matrix(x_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "008890e6-0978-406d-a945-f392a85452fb",
   "metadata": {},
   "source": [
    "---\n",
    "Compare to solution given by least square.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e46db3d-d55a-4676-8ce6-958a1aca756d",
   "metadata": {},
   "outputs": [],
   "source": [
    "la.lstsq(A,b)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb490963-9694-4d9e-bce1-61872843a84d",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "Compare to pseudoinverse to pseudoinverse given by numpy.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c177571-adf5-4b61-8c20-c7bf499eed8b",
   "metadata": {},
   "source": [
    "## Your turn...\n",
    "Answer the following 4 questions.  Go to gradescope Lecture #10 and enter your answers.\n",
    "\n",
    "1. Upload your rank 5 image of MonaLisa.\n",
    "2. What is the largest singular value for the matrix $A= \\begin{bmatrix}1 & 2 & 1 \\\\ 1 & 1 & -1 \\end{bmatrix}$?\n",
    "3. Use python to solve the system using pseudo inverse. You can use la.pinv.  $$ \\begin{bmatrix}1 & 2 & 1 \\\\1 & 1 & -1 \\end{bmatrix}\\ \\begin{bmatrix}x_1  \\\\ x_2\\\\x_3  \\end{bmatrix}= \\begin{bmatrix}1  \\\\ 1  \\end{bmatrix} $$ \n",
    "4. The image of Bernhard Reimann is (246,225) pixels.  So there are 246 * 225 = 55,350 values necessary to store or transmit the image.  How many values are necessary for the rank 5 approximation? <br> Hint: U.shape=(246,246) and V.shape = (225,225).  Also the answer is less than 5000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06ff228b-3d14-45de-b144-a52e5272d0f2",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}