{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2b3eb4ce",
   "metadata": {},
   "source": [
    "Quiz 7, MCS 471 on Friday 7 October 2022, due 10:50am."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5837155c",
   "metadata": {},
   "source": [
    "Fit the data points\n",
    "$(0, 0.449)$, $(1, 0.845)$, $(2, 0.333)$, $(3, -3.63)$\n",
    "by a quadratic polynomial $q$ in the least squares sense.\n",
    "Do the following:\n",
    "\n",
    "1. Setup the linear system that must be solved to compute the coefficients of $q$.\n",
    "\n",
    "2. Illustrate how the QR decomposition applies to solve the system.\n",
    "\n",
    "3. Explain the verification of the least squares property of the computed solution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76167171",
   "metadata": {},
   "source": [
    "# Answer to Question 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ded9da40",
   "metadata": {},
   "source": [
    "The right hand side vector of the linear system are the second coordinates of the data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9caf0985",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Vector{Float64}:\n",
       "  0.449\n",
       "  0.845\n",
       "  0.333\n",
       " -3.63"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = [0.449, 0.845, 0.333, -3.63]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20cc909c",
   "metadata": {},
   "source": [
    "We store the first coordinate of each point in ``x``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "339d1540",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Vector{Float64}:\n",
       " 0.0\n",
       " 1.0\n",
       " 2.0\n",
       " 3.0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = [0.0, 1.0, 2.0, 3.0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "941c7782",
   "metadata": {},
   "source": [
    "With ``x`` we define the Vandermonde matrix.  The $i$th row of this matrix is\n",
    "\n",
    "$$\n",
    "    1, x_i, x_i^2\n",
    "$$\n",
    "\n",
    "and we stop at exponent 2 because we want to fit with a quadric polynomial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "72357f92",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×3 Matrix{Float64}:\n",
       " 1.0  0.0  0.0\n",
       " 1.0  1.0  1.0\n",
       " 1.0  2.0  4.0\n",
       " 1.0  3.0  9.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V = zeros(4, 3)\n",
    "V[:,1] = ones(4,1)\n",
    "V[:,2] = x\n",
    "V[:,3] = x.^2\n",
    "V"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b20da098",
   "metadata": {},
   "source": [
    "The ``.`` before the ``2`` in ``x.^2`` indicates a componentwise operation, that is: we square every element of ``x`` instead of computing ``x*x``."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99c53164",
   "metadata": {},
   "source": [
    "# Answer to Question 2 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7dea667d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}, Matrix{Float64}}\n",
       "Q factor:\n",
       "4×4 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}:\n",
       " -0.5   0.67082    0.5   0.223607\n",
       " -0.5   0.223607  -0.5  -0.67082\n",
       " -0.5  -0.223607  -0.5   0.67082\n",
       " -0.5  -0.67082    0.5  -0.223607\n",
       "R factor:\n",
       "3×3 Matrix{Float64}:\n",
       " -2.0  -3.0      -7.0\n",
       "  0.0  -2.23607  -6.7082\n",
       "  0.0   0.0       2.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "Q, R = qr(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ef890d7",
   "metadata": {},
   "source": [
    "We apply the QR factorization to the system $V c = y$, where $c$ are the coefficients of the quadratic polynomial in the data fit."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64745f3c",
   "metadata": {},
   "source": [
    "After the QR factorization, the system $V c = y$ becomes $QR c = y$ and we then apply the orthogonality property of $Q$ to compute the right hand side vector\n",
    "in the system\n",
    "\n",
    "$$\n",
    "    R c = Q^T y.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ad40a27",
   "metadata": {},
   "source": [
    "One difficulty is that the ``Q`` in the output of ``qr`` is a 4-by-4 matrix and we want to work with a 4-by-3 matrix.  So we work with the first three columns of ``Q``, defining the matrix ``QQ`` first, before we compute the right hand side vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ea65ff53",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×3 Matrix{Float64}:\n",
       " -0.5   0.67082    0.5\n",
       " -0.5   0.223607  -0.5\n",
       " -0.5  -0.223607  -0.5\n",
       " -0.5  -0.67082    0.5"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "QQ = Q[:, 1:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f6ea8e98",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       "  1.0015\n",
       "  2.850763064514482\n",
       " -2.1795"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = QQ'*y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "70adc868",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       "  0.3218499999999995\n",
       "  1.9943499999999994\n",
       " -1.0897499999999998"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = R\\b"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fada0e18",
   "metadata": {},
   "source": [
    "# Answer to Question 3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f23d5586",
   "metadata": {},
   "source": [
    "The residual is perpendicular to the column space of the matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d876520b",
   "metadata": {},
   "source": [
    "We solved the system $V c = y$, we first compute the residual vector, a.k.a. the backward error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "21852dd5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Vector{Float64}:\n",
       "  0.12715000000000048\n",
       " -0.3814499999999992\n",
       "  0.3814500000000008\n",
       " -0.12714999999999943"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = y - V*c"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "819b0a51",
   "metadata": {},
   "source": [
    "Computing the inner product of ``r`` with the columns of ``V`` is one instruction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "1f098e48",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×3 adjoint(::Vector{Float64}) with eltype Float64:\n",
       " 2.66454e-15  4.10783e-15  9.10383e-15"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r'*V"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdab9dc2",
   "metadata": {},
   "source": [
    "And we see that all numbers are around the machine precision.  Doing the verification in one instruction happens below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "6c007d46",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0337005117659315e-14"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm((y - V*c)'*V)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27b27a3a",
   "metadata": {},
   "source": [
    "This last number verifies the least squares property of the computed coefficients of the quadratic fit of the data."
   ]
  }
 ],
 "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}