{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 18.06 pset 9\n",
    "\n",
    "Due Wednesday, November 7 at 10:55am"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 1 (15 points)\n",
    "\n",
    "**(a)** If $Q$ is an orthogonal matrix ($Q^T = Q^{-1}$), explain why it follows from the rules for determinants that $\\det Q$ must be ........ or ........?\n",
    "\n",
    "**(b)** If $P$ is a $3\\times 3$ projection matrix onto a 2d subspace, then its determinant must be ........?\n",
    "\n",
    "**(c)** The following code generates 20 random 5×5 \"anti-symmetric\" (or \"skew-symmetric\") matrices, and prints their determinants. This is a square matrix $A$ with $A^T=−A$. Explain what you observe using the properties of determinants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 5  # you can try changing this too if you want\n",
    "\n",
    "for i = 1:20\n",
    "    A = randn(m,m) # a random m×m matrix\n",
    "    A = A - A' # make it skew-symmetric by subtracting its transpose\n",
    "\n",
    "    println(round(det(A), 13)) # print determinant, rounded to 13 digits\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 2 (10 points)\n",
    "\n",
    "The $2\\times 2$ matrix $A = \\begin{pmatrix} 1 & 4 \\\\ 2 & 3 \\end{pmatrix}$ has eigenvalues $\\lambda_1 = 5$ and $\\lambda_2 = -1$, with corresponding eigenvectors $x_1 = (1,1)$ and $x_2 = (-2,1)$.\n",
    "\n",
    "Find the eigenvalues and eigenvectors of $B = 2A + 3I$.\n",
    "\n",
    "(Before you jump into solving quadratic equations, think about what happens if you multiply $B$ by $x_1$ or $x_2$.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 3 (10 points)\n",
    "\n",
    "(Based on Strang, section 6.2, problem 5.)\n",
    "\n",
    "**(a)** If the eigenvectors of $A$ are the columns of $I$ then A is a ........ matrix.\n",
    "\n",
    "**(b)** If the eigenvector matrix $X$ is upper triangular, then why must $A$ also be upper triangular?  (Note: the inverse of an upper-triangular matrix is upper triangular.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 4 (10 points)\n",
    "\n",
    "**(a)** Show that the trace of $A^T A$ must always be $\\ge 0$ by deriving a simple formula for $\\mbox{trace}(A^T A)$ in terms of the matrix entries $a_{ij}$ (i-th row, j-th column) of $A$.  This is called the *Frobenius norm* $$\\Vert A \\Vert_F = \\sqrt{\\mbox{trace}(A^T A)}$$ of the matrix.\n",
    "\n",
    "Check this in Julia for a simple matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = [1 2\n",
    "     3 4]\n",
    "trace(A'A) # trace of AᵀA should match your simple formula"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b)** Using the SVD $A = U\\Sigma V^T$, derive a simple relationship between the Frobenius norm $\\Vert A \\Vert_F$ and the singular values $\\sigma_1, \\ldots, \\sigma_r$ of $A$.  (The identity $\\mbox{trace}(BC) = \\mbox{trace}(CB)$ from class will be helpful.)\n",
    "\n",
    "Check this in Julia by computing the singular values of $A$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "svdvals(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 5 (7×4=28 points)\n",
    "\n",
    "(Based on Strang, section 6.2, problem 9.)\n",
    "\n",
    "Suppose we form a sequence of numbers $g_0,g_1,g_2,g_3$ by the rule\n",
    "\n",
    "$$\n",
    "g_{k+2} = (1-w) g_{k+1} + w g_k\n",
    "$$\n",
    "\n",
    "for some scalar $w$.  If $0 < w < 1$, then $g_{k+2}$ could be thought of as a **weighted average** of the previous two values in the sequence.  For example, for $w = 0.5$ (equal weights) this produces the sequence\n",
    "$$\n",
    "g_0,g_1,g_2,g_3,\\ldots = 0, 1, \\frac{1}{2}, \\frac{3}{4}, \\frac{5}{8}, \\frac{11}{16}, \\frac{21}{32}, \\frac{43}{64}, \\frac{85}{128}, \\frac{171}{256}, \\frac{341}{512}, \\frac{683}{1024}, \\frac{1365}{2048}, \\frac{2731}{4096}, \\frac{5461}{8192}, \\frac{10923}{16384}, \\frac{21845}{32768}, \\ldots\n",
    "$$\n",
    "\n",
    "**(a)** If we define $x_k = \\begin{pmatrix} g_{k+1} \\\\ g_k \\end{pmatrix}$, then write the rule for the sequence in matrix form: $x_{k+1} = A x_k$.  What is $A$?\n",
    "\n",
    "**(b)** Find the eigenvalues and eigenvectors of A (your answers should be a function of $w$).  Check your answers with the `λ, X = eig(A)` function in Julia for $w=0.1$.\n",
    "\n",
    "**(c)** What happens to the eigenvalues and eigenvectors as $w$ gets closer and closer to $-1$?  Is there a still a basis of eigenvectors and a diagonalization of $A$ for $w=-1$?\n",
    "\n",
    "**(d)** The eigenvalues immediately tell which of these three possibilities occurs for $0 < w < 1$: the sequence *diverges*, *decays*, or *goes to a nonzero constant* as $n\\to\\infty$?    Does this behavior depend on the starting vector $x_0$?\n",
    "\n",
    "**(e)** Find the limit as $n\\to\\infty$ of $A^n$ (for $0 < w < 1$) from the diagonalization of $A$.  (Your answer should be a function of $w$.  Google the formula for the inverse of a $2\\times 2$ matrix if you need it.)\n",
    "\n",
    "**(f)** For $w=0.5$, if $g_0 = 0$ and $g_1 = 1$, i.e. $x_0 = \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix}$, then show that the sequence approaches 2/3.\n",
    "\n",
    "**(g)** With $w=0.5$, $g_0 = 0$, and $g_1 = 1$ as in the previous part, how fast does $g_n - 2/3$ go to zero?  In particular, you should find that $\\frac{g_{n+1} - 2/3}{g_n - 2/3}$ decays proportional to $\\alpha^n$ for some $\\alpha$.   Check your answer by the using the following Julia code, which numerically computes the sequence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function gsequence(n, w=0.5)\n",
    "    g = [0.0, 1.0]\n",
    "    for i = 3:n\n",
    "        push!(g, (1-w)*g[end] + w *g[end-1])\n",
    "    end\n",
    "    return g\n",
    "end\n",
    "gsequence(25)  # compute gₙ for n=0,1,…,24"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gsequence(25) .- 2/3  # compute gₙ - 2/3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gsequence(25, -1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.2",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}