{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "21a9eae2",
   "metadata": {},
   "source": [
    "# 18.065 Problem Set 3\n",
    "\n",
    "Due Friday, March 17 at 1pm."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a19ccce3",
   "metadata": {},
   "source": [
    "## Problem 1 (4+4+4+4 points)\n",
    "\n",
    "Suppose $x = A^{-1} b$ for an invertible $m \\times m$ square matrix $A$.   If we perturb the right hand-side $b$ to $b + \\Delta b$, the new solution is $x+\\Delta x$ where $\\Delta x = A^{-1} \\Delta b$.\n",
    "\n",
    "In class, I claimed that\n",
    "$$\n",
    "\\frac{\\Vert \\Delta x \\Vert_2}{\\Vert x \\Vert_2} \\le \\kappa(A) \\frac{\\Vert \\Delta b \\Vert_2}{\\Vert b \\Vert_2}\n",
    "$$\n",
    "where $\\kappa{A} = \\sigma_1 / \\sigma_m$ is the condition number of $A$.  In this problem, you will prove that fact.\n",
    "\n",
    "**(a)** Show that $\\Vert \\Delta x \\Vert_2 \\le \\Vert A^{-1} \\Vert_2 \\, \\Vert \\Delta b \\Vert_2$.  (Hint: this should follow easily from the definition of the induced norm, which you studied in pset 2.)   Using the same proof, also show that $\\Vert b \\Vert_2 \\le \\Vert A \\Vert_2 \\, \\Vert x \\Vert_2$\n",
    "\n",
    "**(b)** Why does (a) tell you that $\\frac{\\Vert \\Delta x \\Vert_2}{\\Vert x \\Vert_2} \\le \\Vert A \\Vert_2 \\, \\Vert A^{-1} \\Vert_2 \\frac{\\Vert \\Delta b \\Vert_2}{\\Vert b \\Vert_2}$?\n",
    "\n",
    "**(c)** Why does $\\Vert A \\Vert_2 \\, \\Vert A^{-1} \\Vert_2 = \\sigma_1 / \\sigma_m = \\kappa(A)$, the ratio of the largest to the smallest singular values of $A$?   Why does $\\kappa(A) = \\kappa(A^{-1})$?\n",
    "\n",
    "**(d)** The following code generates a random 1000x1000 matrix with a condition number of $10^6$.   How close does $\\frac{\\Vert \\Delta x \\Vert_2 \\, \\Vert b \\Vert_2}{\\Vert x \\Vert_2\\, \\Vert \\Delta b \\Vert_2}$ get to the upper bound above if you try a random $\\Delta b$ vector for a random $b$?   Change the `????` lines in the code below to choose a *particular* $b$ and $\\Delta b$ such that $\\frac{\\Vert \\Delta x \\Vert_2 \\, \\Vert b \\Vert_2}{\\Vert x \\Vert_2\\, \\Vert \\Delta b \\Vert_2}$ reaches the upper bound."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0c203c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "\n",
    "m = 1000\n",
    "U = qr(randn(1000,1000)).Q # random unitary U\n",
    "V = qr(randn(1000,1000)).Q # random unitary V\n",
    "σ = exp10.(range(0, -6, length=m)) # log-spaced singular vals from 1 to 10⁻⁶\n",
    "A = U * Diagonal(σ) * V' # construct A from its SVD\n",
    "@show cond(A) # condition number ≈ 10⁶\n",
    "\n",
    "b = randn(m) # a random b\n",
    "x = A \\ b\n",
    "\n",
    "Δb = randn(m) * 1e-3 # random perturbation\n",
    "Δx = A \\ Δb # change in solution\n",
    "\n",
    "(norm(Δx) * norm(b)) / (norm(x) * norm(Δb))  # ‖Δx‖₂‖b‖₂ / ‖x‖₂‖Δb‖₂"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3631ff5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# change this code so that ‖Δx‖₂‖b‖₂ / ‖x‖₂‖Δb‖₂ attains its upper bound.\n",
    "b = ????\n",
    "x = A \\ b\n",
    "Δb = ????\n",
    "Δx = A \\ Δb\n",
    "(norm(Δx) * norm(b)) / (norm(x) * norm(Δb))   # ‖Δx‖₂‖b‖₂ / ‖x‖₂‖Δb‖₂"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b094adc3",
   "metadata": {},
   "source": [
    "## Problem 2 (10 points)\n",
    "\n",
    "In class, I said that solving least-square problems via $\\hat{R}\\hat{x} = \\hat{Q}^* b$ using the thin QR factorization ($A = \\hat{Q} \\hat{R}$ where both $A$ and $\\hat{Q}$ are $m \\times n$ with rank $n$) is good because it avoids squaring the condition number.\n",
    "\n",
    "In particular, I said $\\kappa(A^T A) = \\kappa(A)^2$ but $\\kappa(\\hat{R}) = \\kappa(A)$.  Derive this. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bffaf19f",
   "metadata": {},
   "source": [
    "## Problem 3 (10 points)\n",
    "\n",
    "Another form of regularized least-squares problem is to assume you know something about what the solution $x$ should look like, i.e. you have a \"prior\" $x_p$ guess of the solution.   Suppose you then want to solve\n",
    "$$\n",
    "\\min_x \\left(\\Vert b - Ax \\Vert_2^2 + \\delta^2 \\Vert x - x_p \\Vert_2^2 \\right)\n",
    "$$\n",
    "where $\\delta \\ne 0$ is a regularization strength (larger δ means that the solution will be closer to $x_p$).\n",
    "\n",
    "Show how this is equivalent to solving an equation of the form $(A^T A + ?)\\hat{x} = ?$.  (Very similar to how we handled the Tikhonov regularization in class, which corresponds to the special case $x_p = 0$!)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d59e8a0",
   "metadata": {},
   "source": [
    "## Problem 4 (10 points)\n",
    "\n",
    "In pset 1, problem 3, you took a $1000 \\times 1000$ matrix $A$ of rank 4 and multiplied it by a random $1000 \\times 10$ matrix $X$, and we argued that the rank of $A$ almost certainly equalled the rank of $AX$ or $(AX)^T (AX)$ or $X^T AX$.   However, the *values* of the singular values were not the same.\n",
    "\n",
    "In this, problem, you'll see what happens if you orthonormalize the outputs $AX$.   We can orthonormalize $AX$ using `Q = Matrix(qr(A*X).Q)` in Julia, which forms the thin QR factorization.\n",
    "\n",
    "For the same rank-4 `A` and random `X` as in pset 1, problem 3 (or at least, generated randomly by the same code, reproduced below), compute `Q` and use it to obtain a *much smaller matrix* (e.g. 10x10, 1000x10, or 10x1000) whose 4 biggest singular values *nearly match* those of `A`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74797236",
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "A = randn(1000, 4) * randn(4, 1000) # random 1000x1000 matrix of rank 4\n",
    "@show svdvals(A)[1:6]\n",
    "\n",
    "X = randn(1000, 10)\n",
    "Q = Matrix(qr(A*X).Q)\n",
    "\n",
    "S = ???? # much smaller matrix using Q and/or A and/or X\n",
    "@show svdvals(S)[1:6]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9415ba3",
   "metadata": {},
   "source": [
    "## Problem 5 (10 points)\n",
    "\n",
    "Suppose that $X$ is a random $m \\times n$ matrix (drawn from some probability distribution) and its mean (expectation value) is $E[X] = X_0$.\n",
    "\n",
    "Derive the following identity for the variance (which we used in class):\n",
    "\n",
    "$$\n",
    "E[\\Vert X - X_0 \\Vert_F^2] = E[\\Vert X \\Vert_F^2] - \\Vert X_0 \\Vert_F^2\n",
    "$$\n",
    "\n",
    "(If you just write out the Frobenius norm definition $\\Vert A \\Vert_F^2 = \\text{tr}[A^T A]$, the problem should be very straightforward.  Remember that the expectation value is conceptually just a sum, so it can be moved inside any linear operation like the trace.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55b01fd5",
   "metadata": {},
   "source": [
    "## Problem 6 (5+5 points)\n",
    "\n",
    "Suppose we are solving a generalized/weighted least-square problem:\n",
    "$$\n",
    "\\min_{x \\in \\mathbb{R}^n} \\Vert b - Ax \\Vert_W\n",
    "$$\n",
    "where $A$ is $m \\times n$, $W = W^T$ is an $m \\times m$ positive-definite \"weight\" matrix, and $\\Vert y \\Vert_W = \\sqrt{y^T W y}$ is a weighted $\\ell^2$ norm.\n",
    "\n",
    "**(a)** Show that this is equivalent to an ordinary least-square problem of minimizing $\\Vert d - Cx \\Vert_2 = \\Vert b - Ax \\Vert_W$ for some matrix $C$ and vector $d$.   (Hint: use the fact that any positive-definite matrix $W$ can be factored as $W = $ \\_\\_\\_\\_\\_.)\n",
    "\n",
    "**(b)** Show that the normal equations $C^T C \\hat{x} = C^T d$ are equivalent to the generalized normal equations $A^T W A \\hat{x} = A^T W b$ given in class."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88ab5b5e",
   "metadata": {},
   "source": [
    "## Problem 7 (5+5 points)\n",
    "\n",
    "The most common form of least-squares is linear regression, i.e. fitting $m$ data points $(a_i, b_i)$ to a model of the form $b(a) = x_1 + a x_2$.\n",
    "\n",
    "Suppose the data points $b_i$ have independent random errors with equal variances $\\sigma^2$ (i.e. they are \"homoscedastic\").   (This is the case in which Gauss–Markov says that ordinary least-squares minimizes the variance.)   In this case, many sources give simple explicit formulas for the variances of the fit coefficients $x_1$ and $x_2$\n",
    "\n",
    "$$\n",
    "\\text{variance of }\\hat{x}_1 = \\frac{\\sigma^2 \\sum_{i=1}^m a_i^2} {m \\sum_{j=1}^m (a_j - \\text{mean } a)^2}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\text{variance of }\\hat{x}_2 = \\frac{\\sigma^2 } {\\sum_{j=1}^m (a_j - \\text{mean } a)^2}\n",
    "$$\n",
    "\n",
    "**Derive these formulas** from the more general equation $W = (A^T V^{-1} A)^{-1}$ for the covariance matrix of $\\hat{x}$ that we found in class for weighted least squares (where $V$ is the covariance matrix of $b$).\n",
    "\n",
    "(The formula for the inverse of a 2x2 matrix will be helpful; google it.)"
   ]
  }
 ],
 "metadata": {
  "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
}
