{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0f2a9965",
   "metadata": {},
   "source": [
    "# 18.085 Problem Set 2 Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9e8ae37",
   "metadata": {},
   "source": [
    "## Problem 1 (5+5+5 points)\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$.\n",
    "\n",
    "### Solution:\n",
    "\n",
    "**(a)** Let's try it in Julia first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "feb2f4ea",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n = 4, eigvals = ComplexF64[-5.551115123125783e-17 - 0.5996786059514322im, -5.551115123125783e-17 + 0.5996786059514322im, -4.683753385137379e-17 - 3.6997879530988254im, -4.683753385137379e-17 + 3.6997879530988254im]\n",
      "n = 5, eigvals = ComplexF64[-2.5589944212612695e-16 + 0.0im, 2.7755575615628914e-17 - 2.6074936578703247im, 2.7755575615628914e-17 + 2.6074936578703247im, 1.1102230246251565e-16 - 4.298750473140713im, 1.1102230246251565e-16 + 4.298750473140713im]\n"
     ]
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "\n",
    "for n in (4, 5)\n",
    "    A = randn(n, n)\n",
    "    A = A - A'\n",
    "    println(\"n = $n, eigvals = \", eigvals(A))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3040f46d",
   "metadata": {},
   "source": [
    "Looking closely, we see that the imaginary parts of tiny compared to the real parts — not including roundoff errors, the **eigenvalues are purely imaginary (or zero)**.  Moreover, because $A$ is real, the **eigenvalues come in complex-conjugate pairs** $\\lambda = \\pm i \\alpha$ except for $\\lambda = 0$ (which only appears for **odd n** since the number of eigenvalues equals n).\n",
    "\n",
    "How do we prove this?  Suppose $Ax = \\lambda x$ for a real $A = -A^T =--\\overline{A^T}$ and $x \\ne 0$ (which may be complex!), and follow the hint and consider\n",
    "$$\n",
    "\\overline{x^T} A x = \\lambda \\Vert x \\Vert^2 = \\overline{(A^T x)^T} x = \\overline{(-A x)^T} x = -\\bar{\\lambda} \\Vert x \\Vert^2 \\, ,\n",
    "$$\n",
    "so it follows that $\\boxed{\\lambda = -\\bar{\\lambda}}$, i.e. $\\lambda$ must be purely imaginary (or zero).\n",
    "\n",
    "(The fact that they come in complex-conjugate pairs follows for any real $A$ simply by conjugating both sides of the eigen-equation: $\\overline{Ax = \\lambda x} \\implies A\\bar{x} = \\bar{\\lambda} \\bar{x}$, so $\\bar{\\lambda}$ is also an eigenvalue.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80a0119c",
   "metadata": {},
   "source": [
    "**(b)** We have\n",
    "$$\n",
    "A = \\begin{pmatrix} q_1 & q_2 & q_3 \\end{pmatrix} \\begin{pmatrix} -2 & & \\\\ & 3 & \\\\ & & -1 \\end{pmatrix} \\begin{pmatrix} q_1 & q_2 & q_3 \\end{pmatrix}^T\n",
    "$$\n",
    "which *almost* looks like an SVD except that the singular values must be nonnegative.  So, we just need to flip the signs of the negative eigenvalues, and flip the signs of *one* of the corresponding q's.  There are multiple possible choices of sign:\n",
    "$$\n",
    "A = \\underbrace{\\begin{pmatrix} \\pm q_1 & \\pm q_2 & \\pm q_3 \\end{pmatrix}}_U \\underbrace{\\begin{pmatrix} 2 & & \\\\ & 3 & \\\\ & & 1 \\end{pmatrix}}_\\Sigma \\underbrace{\\begin{pmatrix} \\mp q_1 & \\pm q_2 & \\mp q_3 \\end{pmatrix}^T}_{V^T}\n",
    "$$\n",
    "where the signs on the left match with the corresponding signs on the right.  (For example, either $u_1 = +q_1$ and $v_1 = -q_1$ *or* $u_1 = -q_1$ and $v_1 = +q_1$.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d7bf906",
   "metadata": {},
   "source": [
    "**(c)** Let's see what happens:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "46b9676f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "svdvals(A) = [2.1323056982359665, 1.4445123720679864, 0.9221317940231752]\n",
      "eigvals(B) = [-2.132305698235967, -1.4445123720679862, -0.9221317940231751, -1.0779172922296461e-17, 1.0779172922296461e-17, 0.9221317940231744, 1.4445123720679878, 2.1323056982359665]\n"
     ]
    }
   ],
   "source": [
    "A = randn(5, 3)\n",
    "B = [ 0I A; A' 0I ]\n",
    "@show svdvals(A)\n",
    "@show eigvals(B);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6175ede3",
   "metadata": {},
   "source": [
    "By inspection, the **nonzero eigenvalues of B are ± the singular values of A!**\n",
    "\n",
    "To explain this, we should be able to take a singular value $\\sigma_k$ of $A$ and the corresponding singular vectors $u_k, v_k$, and somehow construct eigenvectors of $B$ for $\\lambda = \\pm \\sigma_k$.   Recall from class that $A v_k = \\sigma_k u_k$, and similarly $A^T u_k = \\sigma_k v_k$.   By inspection of $B$, after a little thought we can see that this means that:\n",
    "$$\n",
    "\\underbrace{\\begin{pmatrix} 0 & A \\\\ A^T & 0 \\end{pmatrix}}_B \\begin{pmatrix} u_k \\\\ \\pm v_k \\end{pmatrix} = \\begin{pmatrix} \\pm A v_k \\\\ A^T u_k \\end{pmatrix} = \\begin{pmatrix} \\pm \\sigma_k u_k \\\\ \\sigma_k v_k \\end{pmatrix} = \\underbrace{\\pm \\sigma_k}_\\lambda \\boxed{ \\begin{pmatrix} u_k \\\\ \\pm v_k \\end{pmatrix} } \\, ,\n",
    "$$\n",
    "so the boxed quantity is, by construction, and eigenvector of $B$ with eigenvalue $\\pm \\sigma_k$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ef45562",
   "metadata": {},
   "source": [
    "What about the zero eigenvalues of $B$?  Those come from $N(A)$ and $N(A^T)$.   In particular, recall that $u_{r+1},\\ldots,u_m$ are an orthonormal basis for $N(A^T)=C(A)^\\perp$ and $v_{r+1},\\ldots,v_n$ are an orthonormal basis for $N(A)=C(A^T)^\\perp$.  For each of these nullpace vectors, we have:\n",
    "$$\n",
    "B \\begin{pmatrix} u_k \\\\ \\vec{0} \\end{pmatrix} = \\begin{pmatrix} \\vec{0} \\\\ A^Tu_k \\end{pmatrix} = 0, \\qquad\n",
    "B \\begin{pmatrix} \\vec{0} \\\\ v_k \\end{pmatrix} = \\begin{pmatrix}  Av_k \\\\ \\vec{0}  \\end{pmatrix} = 0\n",
    "$$\n",
    "i.e. there are $(m-r)+(n-r)=m+n-2r$ orthogonal vectors in $N(B)$.  In this particular case, that corresponds to $(5-3)+(3-3)=2$ vectors in $N(B)$, yielding the two zero (up to roundoff errors) eigenvalues of $B$ seen above from `eigvals(B)`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff2916ae",
   "metadata": {},
   "source": [
    "## Problem 2  (6+6 points)\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": "465c8759",
   "metadata": {},
   "source": [
    "### Solution:\n",
    "\n",
    "**(a)** For the induced norm, we have:\n",
    "$$\n",
    "\\Vert Q_1 A Q_2 \\Vert_2 = \\max_{x\\ne 0} \\frac{\\Vert Q_1 A Q_2 x \\Vert_2}{\\Vert x \\Vert_2} \\\\\n",
    "= \\max_{x\\ne 0} \\frac{\\Vert A \\overbrace{Q_2 x}^y \\Vert_2}{\\Vert x \\Vert_2} = \\max_{y\\ne 0} \\frac{\\Vert A y \\Vert_2}{\\Vert Q_2^T y \\Vert_2}  = \\max_{y\\ne 0} \\frac{\\Vert A y \\Vert_2}{\\Vert y \\Vert_2} = \\Vert A \\Vert_2\n",
    "$$\n",
    "where we have used the fact that $\\Vert Q z \\Vert_2 = \\Vert z \\Vert_2$, from class, for any unitary $Q$ (and note that $Q^T$ is also unitary!), and did a change of variables $y = Q_2 x \\Longleftrightarrow x = Q_2^T y$.\n",
    "\n",
    "For the Frobenius norm, we have:\n",
    "$$\n",
    "\\Vert Q_1 A Q_2 \\Vert_F = \\sqrt{\\text{tr}(Q_2^T A^T Q_1^T Q_1 A Q_2) } \\\\\n",
    "= \\sqrt{\\text{tr}(Q_2^T A^T  A Q_2) } = \\sqrt{\\text{tr}(Q_2 Q_2^T A^T  A) } = \\sqrt{\\text{tr}(A^T  A) } = \\Vert A \\Vert_F \\, ,\n",
    "$$\n",
    "where on the second line we used the cyclic property of the trace to move the $Q_2$ factor to the left and cancel it with $Q_2 Q_2^T = I$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "142f0459",
   "metadata": {},
   "source": [
    "**(b)** From part (a), since U, V are unitary, we have $\\Vert A \\Vert = \\Vert U \\Sigma V^T \\Vert = \\Vert \\Sigma \\Vert$ in both norms ($Q_1 = U$ and $Q_2 = V^T$).   Now, we simply need to calculate them.  Both calculations will involve the $n \\times n$ matrix (for $A$ being $m \\times n$):\n",
    "$$\n",
    "\\Sigma^T \\Sigma = \\begin{pmatrix} \\sigma_1^2 & & & & \\\\\n",
    " & \\ddots & & & \\\\\n",
    " & & \\sigma_r^2 & & \\\\\n",
    " & & & 0 & \\\\\n",
    " & & & & \\ddots \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "For the induced norm, the calculation is very similar to problem 2:\n",
    "$$\n",
    "\\Vert A \\Vert_2 = \\Vert \\Sigma \\Vert_2 = \\max_{x\\ne 0} \\frac{\\Vert \\Sigma x \\Vert_2}{\\Vert x \\Vert_2} = \\sqrt{\\max_{x\\ne 0} \\frac{\\Vert \\Sigma x \\Vert_2^2}{\\Vert x \\Vert_2^2}} = \\sqrt{\\max_{x\\ne 0} \\frac{x^T \\Sigma^T \\Sigma x}{x^T x}} \\\\\n",
    "= \\sqrt{\\max_{x\\ne 0} \\frac{\\sum_{i=1}^r \\sigma_i^2 x^2}{\\sum_{j=1}^n x_j^2} } \\le \\sqrt{\\max_{x\\ne 0} \\frac{\\sum_{i=1}^n \\sigma_1^2 x^2}{\\sum_{j=1}^n x_j^2} } = \\sigma_1\n",
    "$$\n",
    "(a weighted average of $\\sigma_i^2$ bounded above by $\\sigma_1^2$ just as in problem 2)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0d16ac3",
   "metadata": {},
   "source": [
    "For the Frobenius norm case, it's even easier, because we just have the trace of a diagonal matrix $\\Sigma^T \\Sigma$:\n",
    "$$\n",
    "\\Vert A \\Vert_F = \\Vert \\Sigma \\Vert_F = \\sqrt{\\text{tr}(\\Sigma^T \\Sigma) } = \\sqrt{\\sigma_1^2 + \\cdots \\sigma_r^2}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32250444",
   "metadata": {},
   "source": [
    "## Problem 3 (5+5 points)\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.\n",
    "\n",
    "### Solution\n",
    "\n",
    "For both parts, by the Eckhart–Young theorem, we simply need to find the SVD of $A$ and use $\\sigma_1 u_1 v_1^T$ as our rank-1 approximation.  Fortunately, for both of these matrices it is easy to find the SVD by hand.\n",
    "\n",
    "**(a)** Let's find the SVD by looking at eigenvalues and eigenvectors of\n",
    "$$\n",
    "A^T A = \\begin{pmatrix} 4 & 0 \\\\ 0 & 9 \\end{pmatrix} \\, .\n",
    "$$\n",
    "Since this is diagonal, we can read off the eigenvalues (4 and 9) and hence the singular values are $\\sigma_{1} = 3$ and $\\sigma_2 = 2$.   The corresponding eigenvectors are $v_1 = [0,1]$ and $v_2 = [1,0]$.  We also need $u_1 = Av_1 / \\sigma_1 = [3, 0] / 3 = [1, 0]$.  Hence, our best rank-1 approximation is:\n",
    "$$\n",
    "\\sigma_1 u_1 v_1^T = 3 \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix} \\begin{pmatrix} 0 & 1 \\end{pmatrix} = \\boxed{\\begin{pmatrix} 0 & 3 \\\\ 0 & 0 \\end{pmatrix}} \\, .\n",
    "$$\n",
    "\n",
    "**(b)** This is the unitary 2x2 rotation matrix, from class, with $A^T A = I$.   So, its singular values are $\\sigma_1 = \\sigma_2 = 1$, and for its singular vectors we can use *any* unitary matrices with $U = A V \\Sigma = AV$, e.g. $V=I$ and $U=A$ is one possibility.  Since neither singular value is bigger than the other, moreover, there are **infinitely many** choices of \"best rank-1 approximations\" that are all **equally good**.  In particular, we can take *any* unit vector $v_1$ and the corresponding (rotated) unit vector $u_1 = A v_1$.  For example, $v_1 = [1, 0]$ gives:\n",
    "$$\n",
    "\\sigma_1 u_1 v_1^T = \\begin{pmatrix} \\cos\\theta & 0 \\\\ \\sin\\theta & 0 \\end{pmatrix} \\, .\n",
    "$$\n",
    "as one possibility.   (In fact, any single row or column of $A$ works.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d438a8e0",
   "metadata": {},
   "source": [
    "## Problem 4 (4+6+4 points)\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",
    "$$\n",
    "\n",
    "## Solution:\n",
    "\n",
    "**(a)** $x=Qc$ implies $x^T x = (Qc)^T (Qc) = c^T Q^T Q c = c^T c = c_1^2 + \\cdots + c_m^2$.  (We could also do the same thing term-by-term, of course.)   Similarly, $x^T S x = c^T Q^T S Q c$, but $Q^T S Q = \\Lambda$ (eigenvectors **diagonalize** the matrix!), so $x^T S x = c^T \\Lambda c = \\lambda_1 c_1^2 + \\cdots + \\lambda_m c_m^2$.\n",
    "\n",
    "**(b)** From part (a), by expanding an arbitrary $x$ in our eigenvector basis $Q$, we have:\n",
    "$$\n",
    "R(x) = \\frac{\\sum_i \\lambda_i c_i^2}{\\sum_j c_j^2} \\, ,\n",
    "$$\n",
    "which is just a **weighted average** of the λ's!  Intuitively, an average must be smaller than the biggest λ, but we can show this explicitly quite easily.  Since each of the $c_i^2$ is $\\ge 0$, we can make the numerator bigger by replacing each $\\lambda_i$ with $\\lambda_1$, i.e.\n",
    "$$\n",
    "R(x) \\le \\frac{\\sum_i \\boxed{\\lambda_1} c_i^2}{\\sum_j c_j^2} = \n",
    "\\boxed{\\lambda_1} \\frac{\\sum_i  c_i^2}{\\sum_j c_j^2} = \\boxed{\\lambda_1}.\n",
    "$$\n",
    "\n",
    "**(c)** The trick is simply to observe (as we did in class for least-square problems) that maximizing a positive quantity is equivalent to maximizing its square:\n",
    "$$\n",
    "\\Vert A \\Vert_2 = \\max_{x\\ne 0} \\frac{\\Vert A x \\Vert_2}{\\Vert x \\Vert_2} = \\sqrt{\\max_{x\\ne 0} \\frac{\\Vert A x \\Vert_2^2}{\\Vert x \\Vert_2^2} } = \\sqrt{\\max_{x\\ne 0} \\frac{x^T A^T A x}{x^T x} } = \\sqrt{\\max_{x\\ne 0} R(x) }\n",
    "$$\n",
    "The final expression is simply the Rayleigh quotient for $S = A^T A$.   Hence, the maximum $R(x)$ is, by part (b), the biggest eigenvalue of $A^T A$, which is $\\sigma_1^2$.   Then \n",
    "$$\n",
    "\\Vert A \\Vert_2 = \\sqrt{\\sigma_1^2} = \\sigma_1 \\, .\n",
    "$$"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "[encoded data, 284260 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "b945636a",
   "metadata": {},
   "source": [
    "## Problem 5 (6+6 points)\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": 3,
   "id": "722752b9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"/var/folders/q9/2vfl1wt51558xwpl4nv4wt580000gn/T/jl_DIfs34Bc5M\""
      ]
     },
     "execution_count": 3,
     "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": 4,
   "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": 4,
     "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": "markdown",
   "id": "19ea454d",
   "metadata": {},
   "source": [
    "### Solution:\n",
    "\n",
    "**(a)** The main thing is to express the \"prediction\" vector $\\hat{x} = [\\hat{x}_{M+1}, \\ldots \\hat{x}_{N}] = A\\beta$ as a matrix $A$ times the vector $\\beta = [\\beta_1, \\ldots, \\beta_M]$ of coefficients.  By inspection of the formula for $\\hat{x}_{t+1}$ above, this is:\n",
    "$$\n",
    "\\underbrace{\\begin{pmatrix} \\hat{x}_{M+1} \\\\ \\hat{x}_{M+2} \\\\ \\vdots \\\\ \\hat{x}_{N} \\end{pmatrix}}_{\\hat{x}}\n",
    "= \\boxed{\\underbrace{\\begin{pmatrix} x_M & x_{M-1} & \\cdots & x_{1} \\\\\n",
    "x_{M+1} & x_{M} & \\cdots & x_{2} \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots \\\\\n",
    "x_{N-1} & x_{N-2} & \\cdots & x_{N-M} \n",
    "\\end{pmatrix}}_{A}} \\underbrace{\\begin{pmatrix} \\beta_{1} \\\\ \\beta_2 \\\\ \\vdots \\\\ \\beta_M \\end{pmatrix}}_\\beta \\, \n",
    "$$\n",
    "where $A$ is an $\\boxed{(N-M) \\times M}$ matrix.  The \"right-hand-side\" vector $b$ is then simply the actual $x_{t}$ values to compare to the predictions $\\hat{x}$: $\\boxed{b = [x_{M+1}, x_{M+2}, \\ldots, x_N]}$.\n",
    "\n",
    "**(b)** We fill in the code below to construct our $A$ and $b$ from part (a):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e12d9551",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8-element Vector{Float64}:\n",
       "  0.6478100500722224\n",
       " -0.08569044106766723\n",
       "  0.08122941854537981\n",
       " -0.041598056840006246\n",
       "  0.0030540247873219683\n",
       "  0.07540381091515796\n",
       "  0.7940857504360963\n",
       " -0.4748195773198859"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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, N-M, M) # allocate an (N-M)×M matrix A\n",
    "    for i = 1:size(A,1) # fill in A row-by row\n",
    "        t = M-1+i\n",
    "        A[i,:] = x[t:-1:t-M+1] # x[t], x[t-1], …, x[t-M+1]\n",
    "    end\n",
    "    b = x[M+1:N] # x[M+1], x[M+2], …, x[N]\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": "markdown",
   "id": "80fa5517",
   "metadata": {},
   "source": [
    "Now, to get the predictions for the \"test data\" (the predictions for the last 100 points), we can use the same `autoregressive_model` code!  Remember that the predictions are simply $A\\beta$, so we just construct a new \"A\" matrix from the test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e066bd80",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 101876 characters]",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject Text(24.0, 0.5, 'power consumption (MWh)')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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",
    "\n",
    "A_test, b_test = autoregressive_model(test_data, M)\n",
    "predictions = A_test * β\n",
    "\n",
    "plot(predictions, \"ro-\")\n",
    "plot(test_data[end-99:end], \"b*-\")\n",
    "legend([\"predictions\", \"test data\"])\n",
    "title(\"problem 5b validation\")\n",
    "xlabel(\"days\")\n",
    "ylabel(\"power consumption (MWh)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a648bba",
   "metadata": {},
   "source": [
    "The predictions aren't too bad looking!\n",
    "\n",
    "(Note that it is *always* a good idea to add some proper labels to the plot.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88ac965b",
   "metadata": {},
   "source": [
    "## Problem 6 (8 points)\n",
    "\n",
    "Show that the pseudo-inverse $A^+$ obeys $AA^+A=A$ and $A^+ A A^+ = A^+$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76066214",
   "metadata": {},
   "source": [
    "### Solution:\n",
    "\n",
    "We have $A = U\\Sigma V^T$ and $A^+ = V\\Sigma^+ U^T$, so:\n",
    "$$\n",
    "AA^+A = U\\Sigma V^T V\\Sigma^+ U^T U\\Sigma V^T = U\\Sigma \\Sigma^+ \\Sigma V^T\n",
    "$$\n",
    "where the unitary terms cancel.  All we need to do now is to show that $\\Sigma \\Sigma^+ \\Sigma = \\Sigma$ in order to get $A$.  Suppose $A$ is $m \\times n$ with rank $r$.  It might be square, wide, or tall, and it might be rank-deficient.  In general, we have:\n",
    "$$\n",
    "\\Sigma \\Sigma^+ \\Sigma = \n",
    "\\underbrace{\\begin{pmatrix} \\sigma_1 & & & & 0's\\\\\n",
    "                & \\sigma_2 & & & 0's \\\\\n",
    "                & & \\ddots & \\\\\n",
    "                & & & \\sigma_r & 0's\\\\\n",
    "                0's & 0's& 0's & 0's & 0's\\end{pmatrix}}_{\\Sigma: \\; m \\times n}\n",
    "\\underbrace{\\begin{pmatrix} 1/\\sigma_1 & & & & 0's\\\\\n",
    "                & 1/\\sigma_2 & & & 0's \\\\\n",
    "                & & \\ddots & \\\\\n",
    "                & & & 1/\\sigma_r & 0's\\\\\n",
    "                0's & 0's& 0's & 0's & 0's\\end{pmatrix}}_{\\Sigma^+: \\; n \\times m} \\Sigma \\\\\n",
    "= \\underbrace{\\begin{pmatrix} 1 & & & & 0's\\\\\n",
    "                & 1 & & & 0's \\\\\n",
    "                & & \\ddots & \\\\\n",
    "                & & & 1 & 0's\\\\\n",
    "                0's & 0's& 0's & 0's & 0's\\end{pmatrix}}_{\\Sigma\\Sigma^+: \\; m \\times m}\n",
    "  \\underbrace{\\begin{pmatrix} \\sigma_1 & & & & 0's\\\\\n",
    "                & \\sigma_2 & & & 0's \\\\\n",
    "                & & \\ddots & \\\\\n",
    "                & & & \\sigma_r & 0's\\\\\n",
    "                0's & 0's& 0's & 0's & 0's\\end{pmatrix}}_{\\Sigma: \\; m \\times n} = \\Sigma\n",
    "$$\n",
    "so we are done!\n",
    "\n",
    "The proof is nearly the same for $A^+ A A^+$, since\n",
    "$$A^+ A A^+ = V\\Sigma^+ U^T U\\Sigma V^T V\\Sigma^+ U^T = V\\Sigma^+ \\Sigma \\Sigma^+ U^T\n",
    "$$\n",
    "so we simply have to show that $\\Sigma^+ \\Sigma \\Sigma^+ = \\Sigma^+$.  The algebra is identical to above except that we swap $m \\leftrightarrow n$ and $\\sigma \\leftrightarrow 1/\\sigma$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c918611f",
   "metadata": {},
   "source": [
    "## Problem 7 (6+6 points)\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.\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": 7,
   "id": "5f33349c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"/var/folders/q9/2vfl1wt51558xwpl4nv4wt580000gn/T/jl_leVCLcRON8\""
      ]
     },
     "execution_count": 7,
     "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": 8,
   "id": "7be7e0f5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "size(X) = (19968, 80)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "plotimg (generic function with 3 methods)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 9,
   "id": "14b8590f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 125228 characters]",
      "text/plain": [
       "Figure(PyObject <Figure size 630x315 with 80 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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?)\n",
    "\n",
    "### Solution (a):\n",
    "\n",
    "Here, the different \"samples\" (images) are in the *columns*, so we need to take the *mean of all the columns (images)* and subtract this \"mean image\" from each *column*.\n",
    "\n",
    "We can do this in a variety of ways, e.g. by using the `mean(X, dims=2)` function from the built-in `Statistics` module to take the average along the 2nd (column) direction, then subtracting this from `X` with `.-` (broadcasting):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "5524c378",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = A1 = X .- mean(X, dims=2);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18237689",
   "metadata": {},
   "source": [
    "Alternatively, we can apply the projection matrix $P = I - \\frac{oo^T}{o^T o}$ as in class, except that now we need to apply this projection *from the right* to average columns instead of rows.\n",
    "\n",
    "Note also that, as in pset 1, it is critically important to parenthesize things carefully so that we don't end up multiplying large matrices.  We want:\n",
    "$$\n",
    "A=XP^T = XP = X - \\frac{(Xo)o^T}{o^T o}\n",
    "$$\n",
    "so that we only do vector–matrix and vector–vector products, never matrix–matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f3e8b0cf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "o = ones(size(X,2))\n",
    "A2 = X - (X*o)*(o' / o'o);\n",
    "\n",
    "# let's check that it is the same A:\n",
    "A1 ≈ A2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ac9faa8",
   "metadata": {},
   "source": [
    "Now, we'll plot the singular values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "776c2100",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 24336 characters]",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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": "8658a181",
   "metadata": {},
   "source": [
    "As desired, a \"lot\" (well, 13%) of the variance is explained by the first couple of singular values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "62b323f1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1330148445724591"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(σ[1:2].^2) / sum(σ.^2) # fraction of variance in σ₁ and σ₂"
   ]
  },
  {
   "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.\n",
    "\n",
    "### Solution (b):\n",
    "\n",
    "We need a basis vector for the *images*, so we need $\\boxed{u_1 \\text{ and } u_2}$ (which have the same length as an image).  Equivalently, since the samples are the columns, we need to diagonalize the covariance $AA^T = U \\Sigma\\Sigma^T U^T$ (*not* $A^T A$), corresponding to the left singular vectors $U$.\n",
    "\n",
    "Since we want the dot product of each *column* of $A$ with $u_k$, this means that we want $\\boxed{c_k = A^T u_k}$.\n",
    "\n",
    "It doesn't really matter much whether we do the projection on the images in $A$ or $X$, since they only differ by a constant vector (the average image $x_a$).  That is, $A^T u_k$ and $X^T u_k$ differ only by a *constant* $x_a^T u_k$ in each component, so the difference between $X$ and $A$ here is only an *overall shift* in our coordinate system $(c_1, c_2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "4f260943",
   "metadata": {},
   "outputs": [],
   "source": [
    "c1 = A' * U[:,1]; # Aᵀu₁\n",
    "c2 = A' * U[:,2]; # Aᵀu₂"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "16f9fdc0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 91532 characters]",
      "text/plain": [
       "Figure(PyObject <Figure size 2000x1500 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5708fd7",
   "metadata": {},
   "source": [
    "As desired, the blue points (the \"snow-people\") are clustered together at the right, mostly separated from the red points (the \"snowflakes\") to the left.\n",
    "\n",
    "In fact, we **mostly just need `c1`** to discriminate between snow-people and snowflakes.   It is instructive to show this by a histogram of the images vs $c_1$ alone:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "166ceacc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 23296 characters]",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject Text(24.000000000000007, 0.5, '# images')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hist(c1[snowperson_indices], bins=10, color=\"blue\", alpha=0.5, label=\"snow-people\")\n",
    "hist(c1[snowflake_indices], bins=10, color=\"red\", alpha=0.5, label=\"snowflakes\")\n",
    "legend(loc=\"upper left\")\n",
    "title(\"problem 7(b) supplement: c1 only\")\n",
    "xlabel(\"c1\")\n",
    "ylabel(\"# images\")"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.8.0",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
