{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0f2a9965",
   "metadata": {},
   "source": [
    "# 18.085 Problem Set 2\n",
    "\n",
    "Due Friday **March 3** at 1pm."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9e8ae37",
   "metadata": {},
   "source": [
    "## Problem 1\n",
    "\n",
    "**(a)** The eigenvalues of a real *anti-Hermitian* matrix $A = -A^T$ must be \\_\\_\\_\\_\\_\\_\\_\\_\\_\\_\\_\\_.   Derive this by considering $\\overline{x^T} A x$ for an eigenvector $Ax = \\lambda x$, and demonstrate it numerically by constructing a random anti-Hermitian matrix `A = randn(n,n); A = A - A'` in Julia and looking at its eigenvalues `using LinearAlgebra; eigvals(A)` for `n=4` and `n=5`.  (You might want to try the numerical experiment first if you aren't sure what the answer is.)\n",
    "\n",
    "**(b)** Suppose $A$ is a $3 \\times 3$ real-symmetric matrix with eigenvalues $\\lambda_1 = -2$, $\\lambda_2 = 3$, and $\\lambda_3 = -1$, with corresponding orthonormal eigenvectors $q_1, q_2, q_3$.  In terms of these quantities, give the (full) SVD of $A$.\n",
    "\n",
    "**(c)** Construct a random $5 \\times 3$ matrix `A = randn(5, 3)` and form a related *real-symmetric* matrix `B = [ 0I A; A' 0I ]`, corresponding to\n",
    "$$\n",
    "B = \\begin{pmatrix} 0 & A \\\\ A^T & 0 \\end{pmatrix} .\n",
    "$$\n",
    "Compare the eigenvalues of $B$ (`eigvals(B)`) to the singular values of $A$ (`svdvals(A)`).  What do you notice?  Explain it by using the SVD $A = U \\Sigma V^T$ to construct eigenvalues and eigenvectors of $B$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8110f592",
   "metadata": {},
   "source": [
    "## Problem 2\n",
    "\n",
    "**(a)** For any $m \\times n$ real matrix $A$ and any real unitary $m \\times m$ matrix $Q_1$ and any real unitary $n \\times n$ matrix $Q_2$, show that $\\Vert A \\Vert$ = $\\Vert Q_1 A Q_2 \\Vert$ for the norms:\n",
    "$$\n",
    "\\Vert A \\Vert_2 = \\max_{x\\ne 0} \\frac{\\Vert A x \\Vert_2}{\\Vert x \\Vert_2}\\, , \\; \\; \\Vert A \\Vert_F = \\sqrt{\\text{tr}(A^T A) } \\, .\n",
    "$$\n",
    "Do *not* use the relationships of these norms to the singular values of $A$, from class; use only the definitions above.  (Hint: a change of variables may be useful for the first norm, and the cyclic property of the trace for the second.)\n",
    "\n",
    "**(b)** Using the full SVD $A = U \\Sigma V^T$ and the unitary invariance from part (a), show that $\\Vert A \\Vert_2 = \\sigma_1$ and $\\Vert A \\Vert_F = \\sqrt{\\sum_k \\sigma_k^2}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32250444",
   "metadata": {},
   "source": [
    "## Problem 3\n",
    "\n",
    "Find a closest-rank-1 matrix (in the Frobenius norm, for example) to:\n",
    "\n",
    "**(a)** $A = \\begin{pmatrix} 0 & 3 \\\\ 2 & 0 \\end{pmatrix}$\n",
    "\n",
    "**(b)** $A = \\begin{pmatrix} \\cos\\theta & -\\sin\\theta \\\\ \\sin\\theta & \\cos\\theta \\end{pmatrix}$ (where $\\theta$ is some real number).\n",
    "\n",
    "You should be able to do your calculations completely by hand (it's not too hard, honest!), but of course you may use Julia to check your answers if you wish."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d438a8e0",
   "metadata": {},
   "source": [
    "## Problem 4\n",
    "\n",
    "For an $m \\times m$ real-symmetric matrix $S = S^T$, we know that we have real eigenvalues $\\lambda_1 \\ge \\lambda_2 \\ge \\cdots \\ge \\lambda_m$ and can find an orthonormal basis of eigenvectors $Q = \\begin{pmatrix} q_1 & \\cdots & q_m \\end{pmatrix}$.   So, we can write any vector $x \\in \\mathbb{R}^m$ as $x = Qc = q_1 c_1 + \\cdots q_m c_m$ for some coefficient $c$.\n",
    "\n",
    "**(a)** Show that $x^T x = c_1^2 + \\cdots + c_m^2$ and $x^T S x = \\lambda_1 c_1^2 + \\cdots + \\lambda_m c_m^2$.\n",
    "\n",
    "**(b)** Show that the **Rayleigh quotient**\n",
    "$$\n",
    "R(x) = \\frac{x^T S x}{x^T x}\n",
    "$$\n",
    "is *maximized* (over *any* possible $x \\ne 0$, not just eigenvectors) by $R(q_1) = \\lambda_1$.\n",
    "\n",
    "**(c)** If $A$ is any $m \\times n$ real matrix, we know that the squared singular values $\\sigma_i^2$ are the nonzero eigenvalues of $A^T A$.  Use this fact, combined with part (b), to give an alternative proof of why\n",
    "$$\n",
    "\\Vert A \\Vert_2 = \\max_{x\\ne 0} \\frac{\\Vert A x \\Vert_2}{\\Vert x \\Vert_2} = \\sigma_1\n",
    "$$"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "[encoded data, 284260 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "b945636a",
   "metadata": {},
   "source": [
    "## Problem 5\n",
    "\n",
    "The following problem is adapted from the [VMLS book](https://web.stanford.edu/~boyd/vmls/), additional [exercise 13.3](https://web.stanford.edu/~boyd/vmls/vmls-additional-exercises.pdf):\n",
    "\n",
    "![image.png](attachment:image.png)\n",
    "\n",
    "**(a)** Describe the **size** and **entries** of the matrix $A$ and the vector $b$ for which $J = \\Vert A\\beta - b\\Vert^2 / (N-M)$.  This allows you to find the coefficients $\\beta$ that minimize $J$, i.e. the auto-regressive predictor that minimizes the mean-square prediction error, simply by solving a standard least-square problem (`β = A \\ b` in Julia).\n",
    "\n",
    "**(b)** Implement an autoregressive model to predict German residential power demand with data from [this tutorial](https://vitalflux.com/autoregressive-ar-models-with-python-examples/#:~:text=Autoregressive%20models%20are%20based%20on,going%20up%20in%20the%20future.) via the code below.\n",
    "\n",
    "* Fill in the `???` in the `autoregressive_model` to form `A, b` according to (a) for given data `x` and `M`.\n",
    "* Apply it for `M=8` to find the prediction coefficients $\\beta$ from \"training\" data consisting of the power-demand data excluding the last 100 days.\n",
    "* Use the last 100 points as \"test\" data: compare the predictions of your fitted model by plotting the predictions for the last 100 days along with the actual data (fill in the ???)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "722752b9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"/var/folders/w_/rf8p8vtd5jb1d8rglz61jsz40000gp/T/jl_8RI3CLu1qv\""
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# download the raw data (CSV file)\n",
    "import Downloads\n",
    "datafile = Downloads.download(\"https://raw.githubusercontent.com/jenfly/opsd/master/opsd_germany_daily.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3efb4edd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4384×5 Matrix{Any}:\n",
       " \"Date\"            \"Consumption\"     \"Wind\"    \"Solar\"     \"Wind+Solar\"\n",
       " \"2006-01-01\"  1069.18               \"\"        \"\"          \"\"\n",
       " \"2006-01-02\"  1380.52               \"\"        \"\"          \"\"\n",
       " \"2006-01-03\"  1442.53               \"\"        \"\"          \"\"\n",
       " \"2006-01-04\"  1457.22               \"\"        \"\"          \"\"\n",
       " \"2006-01-05\"  1477.13               \"\"        \"\"          \"\"\n",
       " \"2006-01-06\"  1403.43               \"\"        \"\"          \"\"\n",
       " \"2006-01-07\"  1300.29               \"\"        \"\"          \"\"\n",
       " \"2006-01-08\"  1207.98               \"\"        \"\"          \"\"\n",
       " \"2006-01-09\"  1529.32               \"\"        \"\"          \"\"\n",
       " \"2006-01-10\"  1576.91               \"\"        \"\"          \"\"\n",
       " \"2006-01-11\"  1577.18               \"\"        \"\"          \"\"\n",
       " \"2006-01-12\"  1553.28               \"\"        \"\"          \"\"\n",
       " ⋮                                                      \n",
       " \"2017-12-20\"  1559.69             90.88      8.799      99.679\n",
       " \"2017-12-21\"  1520.37            259.039     7.313     266.352\n",
       " \"2017-12-22\"  1423.24            228.773    10.065     238.838\n",
       " \"2017-12-23\"  1272.17            748.074     8.45      756.524\n",
       " \"2017-12-24\"  1141.76            812.422     9.949     822.371\n",
       " \"2017-12-25\"  1111.28            587.81     15.765     603.575\n",
       " \"2017-12-26\"  1130.12            717.453    30.923     748.376\n",
       " \"2017-12-27\"  1263.94            394.507    16.53      411.037\n",
       " \"2017-12-28\"  1299.86            506.424    14.162     520.586\n",
       " \"2017-12-29\"  1295.09            584.277    29.854     614.131\n",
       " \"2017-12-30\"  1215.45            721.247     7.467     728.714\n",
       " \"2017-12-31\"  1107.11            721.176    19.98      741.156"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "summary(days) = \"4383-element Vector{Int64}\"\n",
      "summary(consumption) = \"4383-element Vector{Float64}\"\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 144568 characters]",
      "text/plain": [
       "Figure(PyObject <Figure size 800x1200 with 2 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "(1.0, 200.0)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# extract the data and plot it:\n",
    "\n",
    "using DelimitedFiles, Dates\n",
    "data = readdlm(datafile, ',')\n",
    "display(data)\n",
    "dates = parse.(Date, data[2:end, 1])\n",
    "days = Dates.value.(dates .- dates[1])\n",
    "consumption = Float64.(data[2:end, 2])\n",
    "@show summary(days)\n",
    "@show summary(consumption)\n",
    "\n",
    "using PyPlot\n",
    "\n",
    "figure(figsize=(8,12))\n",
    "subplot(2,1,1)\n",
    "title(\"German power consumption data\")\n",
    "plot(days, consumption, \"r.-\")\n",
    "xlabel(\"days (from 2006-01-01)\")\n",
    "ylabel(\"power consumption (MWh)\")\n",
    "\n",
    "subplot(2,1,2)\n",
    "title(\"German power consumption data (first 200 days)\")\n",
    "plot(days, consumption, \"r.-\")\n",
    "xlabel(\"days (from 2006-01-01)\")\n",
    "ylabel(\"power consumption (MWh)\")\n",
    "xlim(1,200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e12d9551",
   "metadata": {},
   "outputs": [],
   "source": [
    "# fill in the ??? to get A and b from x and m\n",
    "function autoregressive_model(x, M)\n",
    "    N = length(x)\n",
    "    A = similar(x, ???, ???) # ??? x ??? matrix of same type as x\n",
    "    for i = 1:size(A,1) # fill in A row-by row\n",
    "        t = M-1+i\n",
    "        A[i,:] = x[???:???] # hint: each row of A comes from a slice of x\n",
    "    end\n",
    "    b = ???\n",
    "    return A, b\n",
    "end\n",
    "\n",
    "training_data = consumption[1:end-100]\n",
    "\n",
    "# construct your A and b from M=8 and find β\n",
    "M = 8\n",
    "A, b = autoregressive_model(training_data, M)\n",
    "β = A \\ b # least-square fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e066bd80",
   "metadata": {},
   "outputs": [],
   "source": [
    "# test data from last 100+M points:\n",
    "test_data = consumption[end-99-M:end]\n",
    "\n",
    "# plot the predictions for the last 100 points along with test data\n",
    "predictions = ???\n",
    "plot(predictions, \"ro-\")\n",
    "plot(test_data[end-99:end], \"b*-\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88ac965b",
   "metadata": {},
   "source": [
    "## Problem 6\n",
    "\n",
    "Show that the pseudo-inverse $A^+$ obeys $AA^+A=A$ and $A^+ A A^+ = A^+$."
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "[encoded data, 1646012 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "c918611f",
   "metadata": {},
   "source": [
    "## Problem 7\n",
    "\n",
    "In this problem you will use **principal components analysis (PCA)** and the SVD to classify a images.\n",
    "\n",
    "The code below downloads a [collection of 80 monochrome/grayscale drawings](https://math.mit.edu/~stevenj/snow.zip) of [snowpeople](https://en.wikipedia.org/wiki/Snowman) and [snowflakes](https://en.wikipedia.org/wiki/Snowflake) found by a google image search in Spring 2023, shown in miniature here:\n",
    "\n",
    "![image.png](attachment:image.png)\n",
    "\n",
    "Each image is resized to a common size of $156 \\times 128$, and then reshaped into a vector of $156\\times 156 = 19968$ numbers (the grayscale intensity of each pixel), and stored in a matrix $X$ of \"measurements\" for 80 images:\n",
    "$$\n",
    "X = 19968 \\times 80 \\mbox{ matrix} = \\mbox{(data per image)} \\times \\mbox{(# images)}\n",
    "$$\n",
    "Think of this as 80 points in a **19968-dimensional space**, i.e. each image is a point in $\\mathbb{R}^{19968}$!!\n",
    "\n",
    "PCA allows us to figure out which combinations of these 19968 variables explain most of the variations, and allows us to *project* from $\\mathbb{R}^{19968}$ to a *much* lower-dimensional space in order to help us *classify* the data and identify *clusters* of similar images.   In this problem (parts a and b, below), you will perform the PCA clustering analysis yourself.\n",
    "\n",
    "You will need to install some Julia packages to help extract the image data.  \n",
    "```jl\n",
    "import Pkg\n",
    "Pkg.add(\"ZipFile\")\n",
    "Pkg.add(\"Images\")\n",
    "Pkg.add(\"FileIO\")\n",
    "```\n",
    "(or type `] add ZipFile Images FileIO`).\n",
    "\n",
    "You can then run the following code to download the data and reshape it into our $X$ matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5f33349c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"/var/folders/q9/2vfl1wt51558xwpl4nv4wt580000gn/T/jl_aH2TS1Y9Mx\""
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using ZipFile, Images, FileIO, Statistics, LinearAlgebra, PyPlot\n",
    "\n",
    "# download the files as a zip archive:\n",
    "imagezip = download(\"https://math.mit.edu/~stevenj/snow.zip\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7be7e0f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "rows, cols = 156,128  # common size to resize all of the images to\n",
    "\n",
    "# extract the images from the zip file, load them, resize them, and stick them into a big array Xdata\n",
    "Xdata = Float64[]\n",
    "r = ZipFile.Reader(imagezip)\n",
    "for f in r.files\n",
    "    # load image and rescale to rows x cols\n",
    "    image = imresize(load(IOBuffer(read(f))), rows, cols)\n",
    "    imagedata = Float64.(channelview(image))\n",
    "    if length(imagedata) == rows*cols # grayscale\n",
    "        # convert to RGB\n",
    "        bw = vec(imagedata)\n",
    "        imagedata = permutedims([bw bw bw])\n",
    "    elseif size(imagedata, 1) == 4 # RGBA\n",
    "        α = reshape(imagedata[4, :, :], 1, rows, cols) # alpha (transparency)\n",
    "        # convert to RGB, with transparent pixels converted to white:\n",
    "        imagedata = imagedata[1:3, :, :] .* α .+ (1 .- α)\n",
    "    end\n",
    "    # convert back to grayscale (average RGB channels) and collect in Xdata\n",
    "    append!(Xdata, reshape(mean(imagedata, dims=1), rows, cols))\n",
    "end\n",
    "\n",
    "# array of file names, and indices of snowflakes vs snowpeople\n",
    "names = [f.name for f in r.files]\n",
    "snowperson_indices = findall(startswith(\"s\"), names)\n",
    "snowflake_indices = findall(startswith(\"f\"), names)\n",
    "\n",
    "# reshape into 19968x80 matrix X\n",
    "X = reshape(Xdata, :, length(r.files)) # reshape to (imagesize) x (numimages) matrix\n",
    "@show size(X)\n",
    "\n",
    "# a useful function to reshape a vector into a BW image and plot it\n",
    "function plotimg(imgdata, rows=rows, cols=cols)\n",
    "    u = copy(imgdata)\n",
    "    u .-= minimum(u)\n",
    "    u ./= maximum(u)\n",
    "    imshow(reshape(u, rows, cols), cmap=\"gray\")\n",
    "    axis(\"off\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14b8590f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# first, let's plot all of the images to make sure they look okay:\n",
    "figure(figsize=(6.3, 3.15))\n",
    "r, c = 7, 14\n",
    "for i = 1:r, j=1:c\n",
    "    k = (i-1)*c + j\n",
    "    k > size(X,2) && break\n",
    "    subplot(r, c, k)\n",
    "    plotimg(X[:,k])\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4917563f",
   "metadata": {},
   "source": [
    "**(a)** Fill in the code below to compute the matrix `A` with the *mean of each pixel* (over all images) subtracted, and take its SVD to plot the singular values.\n",
    "\n",
    "(Be careful: some treatments of PCA put the different samples in the rows of $X$ and $A$, and some put the samples in the columns.  Which are we doing here?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "776c2100",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = ??? # subtract means from X\n",
    "\n",
    "U, σ, V = svd(A);\n",
    "\n",
    "plot(σ, \"b.-\")\n",
    "title(\"Problem 7(a): Singular values of A\");\n",
    "xlabel(L\"k\");\n",
    "legend([L\"\\sigma_k\"]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52307e61",
   "metadata": {},
   "source": [
    "**(b)** Compute the projections `c1[i]` and `c2[i]` of each image\n",
    "onto the dominant two principal components (vectors u or v in this case?).  That is, reduce the 19968-dimensional data to **2 dimensions**.\n",
    "\n",
    "(Does it matter much whether you do the projection to the images in $A$ or $X$?  Why or why not?)\n",
    "\n",
    "Run the code below to plot `c2` vs `c1` for every image, coloring the points corresponding to snowpeople blue and the points corresponding to snowflakes red.  Hopefully, in this 2d projection, you will see\n",
    "different types of images separated into mostly distinct coordinate regions, with similar images clustered together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f260943",
   "metadata": {},
   "outputs": [],
   "source": [
    "c1 = ???\n",
    "c2 = ???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16f9fdc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "figure(figsize=(20,15))\n",
    "plot(c1[snowperson_indices], c2[snowperson_indices], \"b.\")\n",
    "plot(c1[snowflake_indices], c2[snowflake_indices], \"r.\")\n",
    "text.(c1, c2, string.(1:length(c1)), fontsize=12);\n",
    "\n",
    "title(\"Problem 7(b): Compressed data, visualized on a plane\", fontsize=20);\n",
    "xlabel(L\"$c_1$\", fontsize=20);\n",
    "ylabel(L\"$c_2$\", fontsize=20);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.3",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
