QR.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Plots.PyPlotBackend()"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Random\n",
    "using LinearAlgebra\n",
    "using Statistics\n",
    "using Plots\n",
    "using LaTeXStrings\n",
    "pyplot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6×4 Array{Float64,2}:\n",
       " -0.890685   0.764352    -0.244713    0.922554\n",
       "  1.09337    0.00810613   0.754116   -0.871233\n",
       "  1.36687    1.44676      0.582838    1.27395 \n",
       " -0.109862   1.23566      0.981122   -1.00478 \n",
       "  1.1635    -0.883671    -0.0394113   0.463097\n",
       "  0.458859   0.43612      0.693942   -0.874769"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# generate random data matrix\n",
    "n,d = 6,4\n",
    "X = randn(n,d)\n",
    "\n",
    "# optional: give it linearly dependent columns\n",
    "# X[:,3] = X[:,2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Understanding the pseudoinverse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×6 Array{Float64,2}:\n",
       " -1.26168   0.00544988   0.715924  -0.900509  -1.1001    0.158552\n",
       " -1.20029  -0.246117     0.849552  -0.820945  -1.83688   0.187011\n",
       "  3.03502   0.567742    -1.36888    2.44695    3.70982  -0.204841\n",
       "  1.36334   0.047329    -0.327362   0.745252   1.55705  -0.260946"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# form pseudoinverse\n",
    "Xd = pinv(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Array{Float64,2}:\n",
       "  1.0          -1.25141e-15  -1.04455e-15  -4.57249e-16\n",
       " -6.39156e-16   1.0          -4.04297e-16  -2.45006e-16\n",
       "  3.0523e-16    1.42441e-15   1.0           1.36557e-15\n",
       "  6.93224e-16   3.62154e-16   7.95279e-16   1.0        "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# X†X ≈ I_d\n",
    "Xd*X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6×6 Array{Float64,2}:\n",
       "  0.721359   -0.288244    0.0446677   0.263312    0.104445   -0.188887 \n",
       " -0.288244    0.390872    0.0425729   0.204752    0.22337     0.247742 \n",
       "  0.0446677   0.0425729   0.992797   -0.0430009  -0.0153952   0.0354588\n",
       "  0.263312    0.204752   -0.0430009   0.736462   -0.0736144   0.274883 \n",
       "  0.104445    0.22337    -0.0153952  -0.0736144   0.918078   -0.0935505\n",
       " -0.188887    0.247742    0.0354588   0.274883   -0.0935505   0.240431 "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# XX† !≈ I_n\n",
    "X*Xd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q,R = qr(X)\n",
    "Q = Q[:,1:d];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6×4 Array{Float64,2}:\n",
       " -1.11022e-16   1.11022e-16  -3.05311e-16   1.11022e-16\n",
       "  0.0           1.19696e-16   3.33067e-16  -3.33067e-16\n",
       "  0.0          -2.22045e-16   3.33067e-16  -8.88178e-16\n",
       "  1.38778e-17   0.0           1.11022e-16  -2.22045e-16\n",
       "  0.0          -1.11022e-16   8.32667e-17   5.55112e-17\n",
       "  0.0           5.55112e-17   0.0           1.11022e-16"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X - Q*R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6×4 Array{Float64,2}:\n",
       " -0.382108   -0.361369   -0.272415   -0.608732 \n",
       "  0.46906     0.0267269   0.411939   -0.0211324\n",
       "  0.586391   -0.599405   -0.517968    0.146167 \n",
       " -0.0471311  -0.547338    0.569154   -0.332754 \n",
       "  0.499148    0.421487   -0.0891341  -0.695222 \n",
       "  0.196852   -0.17939     0.394873    0.116512 "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Array{Float64,2}:\n",
       " 2.33098   0.146637   0.859695   0.0921657\n",
       " 0.0      -2.27021   -0.918873  -0.218209 \n",
       " 0.0       0.0        0.911363  -2.22865  \n",
       " 0.0       0.0        0.0       -0.446499 "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Array{Float64,2}:\n",
       "  1.0          -1.53135e-16  -3.56308e-16  -1.06328e-16\n",
       " -1.53135e-16   1.0           2.17932e-16   9.85739e-17\n",
       " -3.56308e-16   2.17932e-16   1.0          -9.20757e-17\n",
       " -1.06328e-16   9.85739e-17  -9.20757e-17   1.0        "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q'*Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# form data from noisy linear model\n",
    "w♮ = randn(d)\n",
    "y = X*w♮ + .1*randn(n);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Array{Float64,1}:\n",
       " -0.10036720947429924\n",
       "  0.1177262734257617 \n",
       "  0.4330665815147097 \n",
       "  0.10108585641128666"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# solve least squares problem to estimate w\n",
    "w = R \\ (Q'*y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.22139501473072573"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# how good is our estimate?\n",
    "norm(w - w♮)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0054774834985621"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute mean square error\n",
    "mean((y - X*w).^2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.9399232672427924e-16"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# let's use the shorthand\n",
    "w_backslash = X \\ y\n",
    "norm(w_backslash - w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Array{Float64,1}:\n",
       " -0.10036720947429938\n",
       "  0.11772627342576154\n",
       "  0.43306658151471006\n",
       "  0.1010858564112867 "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w_backslash"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.2.0",
   "language": "julia",
   "name": "julia-1.2"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}