Householder_QR_NeklyudovSA.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from numpy.testing import assert_allclose"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# QR decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Any real $m \\times n$ matrix can be decomposed into \n",
    "$$\\mathbf{A} = \\mathbf{Q} \\cdot \\mathbf{R}, $$\n",
    "where $\\mathbf{Q}$ is an $m \\times m$ orthogonal matrix ($\\mathbf{Q}^T \\mathbf{Q} = \\mathbf{1}$, here $\\mathbf{1}$ is an identity matrix) and $\\mathbf{R}$ is an $m \\times n$ upper triangular matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $m > n$\n",
    "$$\\mathbf{A} = \\mathbf{Q} \\begin{pmatrix}\n",
    "            \\mathbf{R_1} \\\\\n",
    "            \\mathbf{0}\n",
    "        \\end{pmatrix} = \n",
    "        \\begin{pmatrix}\n",
    "            \\mathbf{Q_1} & \\mathbf{Q_2}\\\\ \n",
    "        \\end{pmatrix}\n",
    "      \\begin{pmatrix}\n",
    "            \\mathbf{R_1} \\\\\n",
    "            \\mathbf{0}\n",
    "       \\end{pmatrix} = \\mathbf{Q_1} \\mathbf{R_1},$$\n",
    "where $\\mathbf{R_1}$ is $m \\times m$ upper triangular, $\\mathbf{Q_1}$ is $m \\times n$ matrix and has orthonormal columns. This operation is called *thin QR* (*reduced*) factorization. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If columns of $\\mathbf{A}$ are linearly independent, then\n",
    "1. The thin factorization $\\mathbf{A} = \\mathbf{Q_1} \\mathbf{R_1}$ is unique;\n",
    "2. Diagonal elements of $\\mathbf{R_1}$ are positive."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way to reduce $\\mathbf{A}$ to $\\mathbf{R}$ and construct $\\mathbf{Q}$ is to implement Householder reflection algorithm column by column."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part I. Construct a Householder reflection of a vector."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given a vector $\\mathbf{x}$, and a plane with a normal vector $\\mathbf{u}$, the Householder transformation reflects $\\mathbf{x}$ about the plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmgAAAG5CAYAAAA6ZcheAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAAFiUAABYlAUlSJPAAAEBbSURBVHhe7d0HmFXVvf7xhSKIoDRB2tBmQFAiPTZsSLNwsV5RgyV2E0sk/rGH2KKXmBtiAXvBgtGgXLGAoAhoTAApAalDG5p0BBxAyp93nb1nzpyzzzAMp6xzzvfzPDyz9xm5lxAjr2ut37sq7N3HAAAAwBmHeF8BAADgCAIaAACAYwhoAAAAjiGgAQAAOIaABgAA4BgCGgAAgGMIaAAAAI4hoAEAADiGgAYAAOAYAhoAAIBjCGgAAACOIaABAAA4hoAGAADgGAIaAACAYwhoAAAAjiGgAQAAOIaABgAA4BgCGgAAgGMIaAAAAI4hoAEAADiGgAYAAOAYAhoAAIBjCGgAAACOIaABAAA4hoAGAADgGAIaAACAYwhoAAAAjiGgAQAAOIaABgAA4BgCGgAAgGMIaAAAAI4hoAEAADiGgAYAAOAYAhoAAIBjCGgAAACOIaABAAA4hoAGAADgGAIaAACAYwhoAAAAjiGgAQAAOIaABgAA4BgCGgAAgGMIaAAAAI4hoAEAADiGgAYAAOAYAhoAAIBjCGgAAACOIaABAAA4hoAGAADgGAIaAACAYwhoAAAAjiGgAQAAOIaABgAA4BgCGgAAgGMIaAAAAI4hoAEAADiGgAYAAOAYAhoAAIBjCGgAAACOIaABAAA4hoAGAADgGAIaAACAYwhoAAAAjiGgAQAAOIaABgAA4BgCGgAAgGMIaAAAAI4hoAEAADiGgAYAAOAYAhoAAIBjCGgAAACOIaABAAA4hoAGAADgGAIaAACAYwhoAAAAjiGgAQAAOIaABgAA4BgCGgA4ZvaiLd4TgGxFQAMAh8xa/KN56LXvzX0vfm82bNnpfQog21TYu4/3DABIoR827DC3/21mUTCrXvUw8/u+eeak42rZdwDZg4AGAA4o3LHb/OavM0zBmkLvk2J9utQ3N/ZuaipVZNMDyBYENABIsd179poHX55jJs/d6H0SrWm9I8w9V7Q0uQ2rep8AyGQENABIsSEjF5sRE1Z6b7FpBU0raVpRA5DZCGgAkEKf/usH85e/L/TeykZn0nQ2TWfUAGQmAhoApMiM/M1mwNDZdovzQNU6spINaZ1b1fQ+AZBJCGgAkAIr1223E5ubt/3sfVI+F53ewFx3XhMGCIAMQ0ADgCTbWrjLhrOgic3y0OCABgg0SAAgM/CvXACQRNrOfOSNeXELZ5K/YpsZ+fWqcm2VAnATK2gAkESD/5FvRn2z2ns7eKe2qW369cihfgPIMAQ0AEgSBTMFtHhom1vdXHNOY9Om2VHeJwAyCQENAJLgu/mb7P2a8diGrFL5UPP2g51MtSoVvU8AZBrOoAFAgum8mc6dHUw4a3D04faH6Fqod79YYZ8BZCYCGgAkkGo0/vDqHDu5WV439W5mXhnQwdzSp5n3ibE3D/iXqgPIPAQ0AEgQO7H5+sFPbB5VtaI59JAK9gYBnT2Tnbv2sIoGZDACGgAkyOD38+1tAQdrzrIt3pMx/XrmeE/GjJy0yhbeAsg8BDQASABtQeqezXiYsbA45GkFzb/eSSt0b40tsM8AMgsBDQDibPLcjeaFj5Z4bwdPW6ThZ9iuO7eJ92TMuKlrzZLVP3lvADIFAQ0A4igeE5tBwrdKVUqrglrR/5/XPl1mnwFkDgIaAMSJpiofePl7W4MRb3OXbvWeQm7s3dQODsjXs9bb654AZA4CGgDESa0jK5n/uamNubtvC3POiceYnLpVvO8cvMhhA3Wind2xjvdmzJCRi70nAJmAmwQAIIG0qjZn6RYzM/9HG7LKu9JVqeIh5oPHTrRffZrgvGHQNFu5IX++tU1RDQeA9EZAA4Ak0vbnwNfm2qufDlRQANPKmSZGRfdy/u9vf2GfAaQ3tjgBIIl0j+bqDdHdZdoSVX2Gvh+LVuIiXda1YdGq2qzFP9rzaADSHwENAJJIW55B5bK6xunxG44z//f4SWboXe3suyY1da7NN3txdEDT9y89q6H3ZsywMfSiAZmALU4ASKIvp60zj785z3sLUW2GQlksqu7Q6pj6zsLv4/SpI63fY1OLutLuuaJliQECAOmHFTQASKK5Ydc2+fZ3sF/ToNoCDQpnUq1KRXPJGQ28N2PeGLMs7j1sAJKLgAYASRR0N+fxzY7ynspP25z+dqi2UHXDAID0RUADgCTRFmRQzUabZkd6T+WnQQENDPhe/nhpUf0GgPRDQAOAJAmawtT2ZfggwME4/5R6tsBWNIww6pvV9hlA+iGgAUCSqKw2UusmB7965tMq2pXdcrw3Y94eu7zEJesA0gcBDQCSJGgFLd7N/5re9FfRNm/7uajEFkB6IaABQBLoPFhQQFP7fzzpAnVdpO77YOIqVtGANERAA4AkmF+wNerQvs6e+atd8aSCW3WricIZ5bVA+iGgAUASqGg2Utu8+K6ehbvu3Cbek7HDAhoaAJA+CGgAkARBAwKtGsdvQCCS7vX0z7dp5Y5VNCC9ENAAIAmSMSAQqV/P4onOT7/9IfAOUABuIqABQIKpnDbyoL6uZ/LPiSWKAuBJx9Wyz7r66a2xrKIB6YKABgAJFnS9Uzz7z0pzTa/G3pOx1z8F3WQAwD0ENABIsKAL0k/ITdyAQDit0p1+Qm37rFU0zqIB6YGABgAJNmNh9IBAvPvPSnPdeU1tP5p8PWt94Hk4AG4hoAFAAv2wYUdUxYWuZGqZU817Szx1remGAZ8uUgfgNgIaACRQrPNnCmnJdFWPxkX/P/VrCvp1AXAHAQ0AEiiooDZZAwLhjqlV2Zx/Sj3vzZgXPlriPQFwEQENABIoKKAla0Ag0mVdGxatounqKZ1HA+AmAhoAJIjOnhWsKfTeiiVzQCCc7v5USPNpolOTnQDcU2HvPt4zACCOtEI18NW53luIai+G3tXOe0s+Feb2e2xqUXHu3X1bmB6d69rnRPjom9Vm3NQ1Jn/lT2b7zt32s/q1DzfnnHiMufzsRvY90vK1headccvNhBnr7c85vNKhpnOrGuahq1t5fwWQ+VhBA4AECbp/M9HXO+2PbjC45IwG3puxtwskahXt3hdmm7/9I9/MXrKlKJzJqvXbzSufLLUhLJLC2X0vfm/GTF5T9HPCfy6QLQhoAJAgQZOSx6doezPcpWc1tNudovs5dU9nvCmcTZm3yXsL9sHEld5TiB/OFOAitWiUvFoSwAUENABIgMIdu82SVT95b8XaNEv+BGckDQpEnkXbuWuP93bwBg1fYMOZtiZvvzjXvHpPB7utGWnjlp9LrKI9+8GiwHAmqZh8BVKJgAYACaDpzcitw5y6VYpWrlJNlRsqsBUNM4z6ZrV9jodWjY80zRtUtWfbeu/7/9OoThV75izIguVb7denRywqWnFTmBt0SxvTt2vojFrNIw8z7fJSuzUMJBtDAgCQAG+MXhZ176UCiw7lu0LnvLTaJTqbNuz+jvZronTv/7X3VExBTsMCjw2bZ98Vzh6/4Tgb6oBsxgoaACSAiwMCkXT9k7+KpqnOERNKngmLN4WxSItWbrMDA6KVMsIZEEJAA4A409Zm0IXkqeo/i0UXqN/Yu6n3pkP7q6LuDY2nvIbRAU107kzn1W69oDnhDPAQ0AAgztTSH3noXmfP/NUql5zaprbtZhOtor37xQr7nAg6mxbLBV3qmzPbHe29ASCgAUCcBdVrtM1za/Us3HXnNvGejB0WSNQqWvsWwVu8nY6tYa47r/jXAICABgBxN3dpaDIxXGmrR6nWuVXNovNxWvl7+ePQmbB40/altjIjqZIkntSnphsM9GP6wuiwDKQDAhoAxFngCppjAwKR+vXM8Z6MGTd1rS2wTQQNAkTasOVn7yk+pi3YbG8w0I/Pp6zxPkWqzPpxmbllxlCz6edt3icoCwIaAMRR/optRfdc+lRd4Z/zcpUCpM6jiYYc3hizzD7H08Ovzw0soo1VTov0NmXTQnPhv58wv/jyDjN0yWjz1/yPvO+gLAhoABBHQdOb6dKC369HyVU0hc140bbpxJnrvbdo2o5EZhi/bpY56+sHTeev7jYfrvqX96kxgxeNMtv3xHe1NJMR0AAgjmYuit7ePCHX3QGBcFrlO6t98SRlZNFueY2fvs58OGmVfQ46gyZLVkdfi4X0ojCmUKZwppAW7oL6J5ovT33EHH5I9BY3ghHQACCOvk+D/rPSXNOrie1Hk69nrbdXVh0MHdh/7sNFZvvO3TacDbmrrfedklauK/SekE527d1tXlv2hd3G1HamtjV9FSscaq5p3NUs6Pac+eCX95h21Zt530FZENAAIE5+2LDD/gini8lb5lTz3tynrjbdMOB77dPyn0VTOLvvxe/tpehyRbdGdpIzaFAgfyUHyNOJtip1pqz1uNvMtdOetoMAPq2S3Znb2yzu/rx5tf1tJq9qfe87OBAENACIk6DVJoUzhbR0cs05TYp+zZpIDZpKDdL73m/tfZvvjFtu6y0UzvwBgNNOqG3v3JTcgCufFOK0FerTHaH3vjCbs2mO0STmEwtGmGZjbjS/m/WKWbgttHUtNQ6rah5oealZ3OMF879tfm0aVQkNnaB8CGgAECezlkQHtHTa3vTVqV7JnH9KPe/NmCEjF3tPsSmQaRtTdLfm3UNmFYUz3cH50NWt7LPkNQxeUfxwYuguUA0U6CL3KfM2cTbNEat3bDIPznnbNPv8JnPv98Psu69e5RrmT8f1sytmj7S+wr7j4BHQACBOZgSUoqbLgECky7o2NFUqhw70a5pT59FKM3VfmApSv/bh5sGrjvXeQnr+sm7gsMDsJVvsCtzwL5bbd90wcNtFze0zUmN54Xq7UqYVs0fnv1eiy0xbl1op04rZPS0usitoiB8CGgDEgbrPCtZEH3RPl4qNSLo79JIzGnhvoVUt9aPFsiLgkL/C2eM3HBd1Abree3Su670F05bon2483ntDsmnrUmfLtGKms2bh9Rhtjmpsz5bNOftpe9aMyczEIKABQBwEndNSbYVKatPVRac3KPr1K3yqG600/qqYhgAUwILCmU8rY327NrIhzqefp1Wz2y/OLbEliuSZvnmxncZsMfZWO52pKU2fpjA1jfmfswbb6UxNaSJxKuzdx3sGAJSTzmmNmBA6Q+VTwLmlT3pXC2i70b+bUxOerwzoUFTD4SoNFuiaJ1FQvLtvC/uM2NRb9tj8983YtTO8T4qdeXQb84djL7NfkTysoAFAHATdIODyBellpZCp7U7R/ZwjvcJZZAaVy548YYAtl40MZyqXnXzGIFswSzhLPgIaABykwh27zfyCrd5bsbZ56TkgEE51GxoY8L37xQqzc9ce7w3pyC+XbT/+Lrud+e3G+d53istltY2p7cxONfK87yDZCGgAcJDmL98adYBe24H+ylO669Olvv3PIxu27IzaykV60EF/XVrul8vqvJnPL5fVwX8NAGgQAKlFQAOAgxRUr5GO/Wex6MzZld2KL1LXKpqmVpEeVI2hSUxVZdwyY2hUuawqMvxyWVr/3cGQAADnqPT03S+W2+t//GuCNOHXuVXNUg986zD76Mk/FP0cFaQ+37+dfU6kAc/PNt/NL9kDpl/n/qok0olWCG8YNK2oSqRfjxxzVU83V1kYEghZt/NH88yiT8zgRaNK9JeJymRvatrTrprRX+YmVtAAOEUhSy30apH3g5boWe3yD78+1/ukJF0LpInD8J+js2GJpuASNCCQrv1nsWgV7brzmnhvxrz/1Uq73Qn3qOVf5bI5o683f5z3bolwpuuX/HLZga36Es4cxgoaAGconPkt8qUZdEsb0y6vuvcWCmcKdJFUdproPi2Fs9v/NtN7C9HZs3cHdvbeMsvNf5lubxaQVNWI3PTUdLMojperf/7Uqd5TetPWpaoy3iz4qkR/mWjr8v6Wl5hf5ZxBf1maYAUNgBO0LeWHs/86tb79Q1OlpUH+7+viMzQKdUHhTBoeHVySGk9Bq2dtmmXW6lm4G89v6j0ZM2rff2c/bNjhvSFVdNj/0smD7OH/WOWyC7o9R7lsmmEFDYAT/HNn2pb8620n2M+Wry001z7xnX0Op/b5N+7raMZPX2ceGzbPfqYW+xt7NzU5davYLVKJXGlLhIGvzo26p1KrSlpdylS/f25W0c0JFMGmjspln8ofaUatnuJ9Uky9ZVox61anrfcJ0g0BDYDTYm1nvXpPB3PXs/+xZ84UzvpflmfObHe0993kufCBf0VNNA69q5295ilTKZwppInOpul2Ab+GA4mnclkd/FdAi6Ry2QEtLjIn1WzpfYJ0xRYnAKflxQg69734fdFAQKrCmSYaI8NZlcqHZnQ4k7a51c2pbWrbZw1JvPzxEvuMxNG25fAVk4rKZcPDmbYtdbZs2pl/sduZhLPMQEAD4LRY1yWtWr/dftWF26kIZ5Lp/WelUc2Gb8LM9UWDA4gvlcu+tPRze77s8ilPRZXL3ty0py2XHdbhTnveDJmDgAbAab1Pqec9RdMQQXj1Q7LNWRY9IHBCbnYENK0SntW+OBi/9tky7wnxoGDml8veMP25qHJZ9ZepKmNI25spl81QBDQAztNQQCQV1/7pxuO9t9TwD8qHy7T+s9Jc06uJPYMm336/IfD3AwdGnWUD5w439T+71naZqdPMd3Slo8wfjr3MLO7+vO0yU9ksMhcBDYDzau0LY5EKd6T2wu61m3dGVUzoYvFsCmgaDDi7Yx3vzZhhowu8JxyoonLZMdHlsgpiCmQFPV+iXDaLENAAOK/WUdGXjm/fudvWcKTKrEU/ek/FWuZUsyEtm1xzTpOi/8xaQZs8d6N9Rtlo61JbmNrK1Jbm1l2hs5WirUtdXK6tTG1p6swZsgcBDYDT3hm33EycWbJnzDdtQeq21GYuyt4BgXB1qlcy54edE3z5k6XeE0qjw/469K/D/xoC0Jkznw77v9f5bnv4X+WyBLPsREAD4CwV0b49NvbVT3MDDuknS9ANAtkyIBDpsq4Nbb2IaJozsrgXxSatn2NrMlSXodqM8NZ/lct+dOL9ti7jkgan0Pqf5QhoAJyk7ctXPllqtzJVRBtkRYq2ONV9FlQrkU3nz8Lp7tFLzii+OeGFj5bYfjQU+2zNNHPW1w+a0ybdZ4tmw6lc9stTH7E/zq/XyfsU2Y6ABsBJfx6+oKjrTFc4NW8QfTB6pff9ZAuaVlTtRLUqFb237KOrrfz//CvXbTfjpq61z9lMq2MKY1otO+efD0eVy/Zt2KWoXFarZ0A4AhoA59z7wmwze0loC/G0E2rbLrSgSU7dJKA7PJNt7tKt3lMxtetnM4WzK7o18t6MeWPMMrNzV2onbVNFwUyXlut8mbYzI8tlda5M58ve6dSfclnEREADkFJ3Pj3TdO//tXn49bl2W1PhbMq8UPeT+s8eurqVfc5rWM1+jRQ+QPDRN6vt3Z0vf5zYg+pBK2ixbjzIJn261LfbnaIKklH7/vvIJn65bM7o6821054uUS6rYOaXy2oyk3JZ7A8BDUBK+duUClrXPvFdUTjTubPHbzjOPkvHY4NLOSfOXGeDnQYKdPZJF6v/J2DCMl60KjS/IGAFLS87BwTCqW4jfBXt3S9WZMUqml8uq2AWWS6rzjKVy67q9SrlsjggBDQAKaPtSf/C83AKZ7oAvVGdKt4nxrTLqx54Dk0/X8HusWHz7ECBVt1+37eF99340/Rm5AF4Fbb6K0fZTpUb+v2QDVt2muHjYk/hpjsFsbtnv26afX6TLZddt7O4G6+oXLYH5bIoHwIagJQJqqrww1nQBeiXn128OhNEAU6rbuHBLt6CCmqzsf8sFl39dGW34ovUP5i4yk69ZhJtXd4280VbLvvnhR+WaP3X1uWL7W4tKpetVjH6mjKgLAhoAFLmp+277Z2aomCmy88fua51YDgTfX77xbnm+KbF57308/Tet2sj83z/dgkNZ+IPL4TL9gGBSLr+qWm9I+yzwtn7X620z+lOh/11tkyH/59Z/ElUueywDnfaw//XN+lOuSwOWoW9+3jPAIBSaGvzwgf+ZQp3FJeLyisDOpicuokNhulGZbUDX51rn3U2bdgDHdN2G/jbjfPNkwtGRPWXSZfarU3/3D62ywyIJ1bQAKCMlqz6KSqcKXQQzqKd2qa27YYTDQpoYCDdjF07w5bLnjxhQFQ461W3vS2WndjlccIZEoKABgBlFFSv0aYZ9Rqx3NKnuONLlRsqsE0Hfrls928GRpXLKoypXPbTkx+iXBYJRUADgDKamR89IHA8AwIx6Wyefz5Pq2hvjS2wzy4qLpf9bVS5rIKZXy6r1n/KZZEMBDQAKKOgqVMGBEp3zTmNvSdjr39ybRXNL5dtMfZWOwAwd2vxVqxfLlvQ8yXKZZF0BDQAKIOCNYW21ytclcqHFp2zQjBVkOg8mmjIQmXCLlA1xqPz37NVGSqXXfLTGu87xeWyqsqgXBapQkADgDKYtZj+s/Lq16O4F03TnfkrinvDkk3lsvd+P8yWyz445+0Srf8KYoOOv9os7v68LZclmCGVCGgAUAZBAS28jw2xaZVR3Wi+F0YlfxVNK2RaKdOK2RMLRkSVyz59wg12xez3eRfQ+g8nENAAoAwCV9Cas4JWVlf1aGxvGZDv5m8KnIhNBJ0p09kynTHTWbPIclmdLdPh/982O5dyWTiFgAYA+6GzZ5GH21W+2roJK2hlpfs5w1fRho1O7ESnpjA1jampTE1nakrTd1LNlnYaU3UZms7UlCbgGgIaAOzHjIXRq2fatlNIQ9ldd16Tot8zraDpPFq8qbdM5bLqMYssl1Vvmcpl/3n6k5TLwnn80wUA9mPuMuo14kG3Lpx/Sj3vzZhhY+K3iqYw1vmru204Cy+XFb9cVuGMclmkCwIaAOxH0Hmp47lBoFwu69rQVKtS0T5rmvPLaevsc3n45bK/+PIOu505ZdNC7zvh5bLPUC6LtERAA4BSbC3cFVgLQcVG+WgV7cLTigtfX/tsqe1HOxA66P/M4k9M63G32QGAWT8u875TXC67oNtzdgCgVbWG3neA9EJAA4BSBN0eoPNn/ioQDtxFpzco+v3T8IVuGCgLVWP8eeGHtirjtpkvmoXbVnnfCZXLPtDy0qJy2aZH1PW+A6QnAhoAlCLo/k1Wzw6OwtkV3Rp5b8a89ulSe1dnLCqTHTh3uC2XvXv261Hlso+0vsKWy+or5bLIFAQ0AChFcEEtAe1g9elS3253ytrNO82ob1bb53DLC9cXlcv+cd67JcpltUKmlTKtmGnljHJZZBoCGgDEoFWd+QVbvbdibfMIaAdLdRvhq2jvfrHCnvcTbV3qbJlWzCLLZXWmTGfLdMZMZ80ol0WmqrB3H+8ZABBGq2e/e+Y/3lvIMbUqmzfv7+S94WBoOODXT35XVALc/fRq5p/1PzWjVk8pUSwrmsLUBeb0lyFbsIIGADEEbW/SfxY/uvrpym7FF6l/+s16M6pgeolw5pfLqseMcIZsQkADgBgYEEgslcs+sn2w2XpEqAut4q7Kptmyk+yzwtjkMwZRLousxRYnAMRw4QP/KjoX5XtlQAeTU7eK94YDpdWx4csnmafyR9r7MqXO+jzTdvYF9rnCIXvNg7+rb05rkGvfgWzFChoABFA5bWQ409Qh4ax8dNB/6JLRtly233d/LQpnsqXOUlO5Vugc2t49Fcz8aXTMAaygAUCAERNWmiEji0OEnNqmthl4bSvvLTsptOr3pqx27tllvtuUb77eMNds3R0KYT/UmWu2HbHBVmNc36S76Z/Xx/ywvIJ58OU55pIzGphLz2xoqlQ+1P61QLYioAFAgMffnBd1T+QtfZrZFvxs9/vnZgXeT1oWG2ouNas6jTM3Nu1uazLC+8s2b/vZVK9KbQYgbHECQIAZC5ngjKVfz+LJywOxq1Kh+e8LqptFPYaaga36RpXLEs6AYgQ0AIigXq4NW3Z6byHacmta/wjvLbspqB5oWK1QYa/5n2s6mAFtzqNcFigDAhoARAjqP1O9hnq7EHKgq2i/6t7YdD62lvcGYH8IaAAQIfj+zSO9J8jG6gXmkNrRv09BOrSsYa7sXr5tUSBbEdAAIELgClpzCmpF1zB1/XiQufOFqWbP+v3/nqiaZMAVLVh9BA4QAQ0AwujsWcGaQu8tRBd7t26SvStoKpd9f+U3pvPIh80Dr8w2h35xii2X3R+Fsvv6tbQhDcCBIaABQJg5S7d4T8VyG1a1IS3bqFz2tWVfmPbvPWr+9MoSU2PC2VHBrGXT2MW915zTmMlXoJwIaAAQJuj+zWwLGQpmf83/yLR59w/mmdfWmHr/6mZqbWzifTekY+ujzNC72plnb+sQ+PvTuVVN07drI+8NwIEioAFAmKAC1uObZcf25qaft5mBc4eblsPvN2++tc00m3Keqbm5+HB/hQrGnNq2hg1mT1z/C7uyKJETnXWqh86dASg/AhoSYu7WFd4TkD4Kd+w2S1b95L0VU8VGJlu9Y5P53axXzAnv/dGM+rsxx07tExXMzu50tHntno5m4FXHFwUzX3gvWujc2bGUzgIHiYCGhHir4Ctz6eRBZuuu0N17QDrQ9ObuPSVvv1MYqVYlMy/vXrhtlbllxlDT7u8Pm4kjjggIZnuLgtk9lx9rGhx9uPedaP4qms6dZXqgBZKBgIaEWL59vZ36Om3SfWZ54XrvU8Bts5dEDwhkYtiYvnmx6ffdX83JI540k0fWML+YfnGJYHbIvj8Zutlg1mm/wcynFbRrejXm3BkQJwQ0JJT+IOj81e/NpPVzvE8Ad81aFFRQmzkBTf87vPDfT5ge/3jOLBjZxLSbfmlUMOvRua55dUBHM6CMwSwcZbRA/BDQkBDbdxffY6jzLWd9/aAd1wdcpa3NoIqNtnnpH9A+WzPN/m/wopGvmFUfH2vazr7AHLm1rvfdksHs7r4tDjiYAYg/AhoSQqEsnIour532tD2IrGfANfMLtpqdu/Z4byHH1Kqc1iWrH676l2k//i5z1ai3TeGYDlHBrOKhxlx0egPz5gOdCWaAYwhoSCp1K/X+9jE7zg+4JKheIx37z/QvQFqtbjH2VnPjp++byl+eFhXMDqtYwQaztx7sbG7p08zWYgBwCwENCVFaANN2y8kTBlDFAafMXbrVeyqWTgMCfrlss89vMv9v7Kem9sTuMYPZmw90ssGMK5gAd1XYu4/3DMSN/pBY8tMa7y1YjcOqmnc69Te96rb3PgFS58IH/mW2Fu7y3kJeGdDB5NSNfZWRC/QvQwpmLyz53OxaWcs0X3pKiVAmhx1awfQ+tb65rGtDQhmQJghoSIiyBDSpWOFQM+j4q82dub29T4Dky1+xzdz8l+neW4iCzLsDO3tv7tE5z6cWjjQvLxlnqqxoYpouO9EcUVjT+24IwQxIXwQ0JESVj/7bbrmU1TWNu5oX291qAxuQbKO+WW0G/yPfews5tU1tM/DaVt6bO1Qu++ziT83zi8aYmqtbBAazqocfanqfUt9ceHp9ghmQpghoSIgKIy/0nsquS+3W5r3O/8/Uq1zD+wRIjsffnGe+nLbOewvRGS2d13LFrB+XmafyR5q3lk00dX44NjCY6caDC0+rb3/dmXr7AZAtCGhIiPIENGlUpbb56MT7TbvqzbxPgMS74uHJZu3m4u4+0YXgkXdOpsKUTQvNY/PfNyNXTjb11xxHMAOyBAENcaernXLGXO+9HbhqFQ83r7a/zVzS4BTvEyBxftiww/zqsSneW0iVyoeaDx490V78nSrj180yf5z3rvlq7fcEMyALUbOBuDvYIlpdsK6L1gfOHe59AiSOLkiPpHqNVIUzv1z27IkDzcKZh5qTp1xrjpvXq0Q4Uxjr1yPHDLu/o7mqZ2PCGZCBCGiIu3jdFKDVAwU1BTYgUWYtiQ5orRpX856Swy+XbT3ut+bibweZDd/XMl3+fZNptbA7wQzIUmxxIu60NaN7/+JF59F0Lk3n04B4+/WT35mCNYXeW8ifb22TlFsENOk8dPFndipz0ZY1ptGqdqZpwYmm0s4jvL8iRJOYqso4/5R6plJF/r0ayAYENMRdvAOaaLJTE56a9ATiRcW0KqgNp63N//vTSQkNQiqXHbpktBmc/5FZU7iFYAYgCv+LR9wlYktSf6D9Y9U/uWgdcRV0/2bLnGoJC0Mql31wztu2yPn+WW+bSoty7VZmy/yzSoQzBTPVfAx7oKMdACCcAdmHFTTEnc7SXDvtae/t4Gjl7I7c3ubmpj3t1VBAPA0ZudiMmLDSewvp27WRue68Jt5bfOhWjcGLRtntzJ279rBiBmC/CGiIu3gENJ0765/bx/Rt1IXbBZAwt/9tppmzdIv3FqLbA3SLQDzM3brCPLlghHmz4Ctjfq64L5i1NY1XdEpYMPvom9Vm3NQ1Jn/lT2b7ztBqc/3ah5tzTjzGXH52I/seafnaQvPOuOVmwoz19uccXulQ07lVDfPQ1e7dogBkEwIa4k4XN/9u1iveW9kpiJ1fr5O5o/n55syj23ifAolRuGO3PX+2e0/JfwSq/+xgpyOnb15sp5BVmVFxV+V9oayj/aHncA2OPtxcdFoDc85Jxxz0itm9L8w2U+Zt8t6i/frcJlEhTeHsvhe/N6vWlzyWcNoJtQloQIqxho6403mx8vjn6U+YD355D+EMSaGVs8hwppsDDiacaUCm+zcDbY/ZqILppvnSU0yXf99ov4aHMwWzu/u2MK8M6GD6dKmf8HAmH0wsuZUbK5xJi0bJrRkBEI2ABmfoDzcgWYIKals3OdJ7OjBaKTt5wgA7vTx+1dwyBbMenevGpQx30PAFNpxpa/L2i3PNq/d0sNuakTZu+dluZfqe/WBRYDiT8v4+AIgfAhri7ocdpf+bfKcaeXYrM5IOUTOliWSZvaTk2TM5oXnZu8/8clmtll347yfMlLVLkxrMfK0aH2maN6hq/+/2PqWeaVSnij1zFmTB8q3269MjFhWtuCnMDbqljR2OkJpHHmba5SW+Aw5A6TiDhrjTgID+4Aqn82UX1D/R9M/rY06q2dIs3LbKtB53W1Qge6dTf9O3YRfvDUgMbW3q/JnOoYV7+6HOpk71St5bMJXLvrT0czM4f5T9+3h/Z8yu7JZjzu5YJ+lXR3Xv/7X3VExBTufQHhs2z74rnD1+w3E21AFwCytoSChVY9yZ29ss6Pacea/z3TacSV7V+oGraE8tHOk9AYkzv2BrVDg7plblUsOZzlZqAKbZmBvNbTNfNMs2bjYtFp2R9BWzslIYi7Ro5TbzyidL7bNWyghngLsIaIg7dT41PaKu+d82vzaLuz9vv+o9kmo0Ik3ZtNBMWj/HewMSI7JaQ2Jd7bRu54/24n6Vy2o6ecOWnbZYVgWzTZZ3LhHMNGSgmo7X7+2YsmDmy9v3awmic2c6r3brBc0JZ4DD2OJE3ClkqcesLP1lnb+62/714bQVqmlOIFEGvjrXfD1rvfcWcsfFubaHzKfWf3WYqVxW25rqLlO5rEpmD9lT8u9tBTNdYh6v/rR4UCfa3/6R772VlIgyXgDxxQoa4k5DAGUtl9WZtEijVk+xZ3uARAm64qmtdzBef+/pHGXO6Ovtluae7YcVrZjpnFl4OPNXzIbe1c6pcCbtWwSvCHY6tgbhDEgDBDSk1CUNTjaNqpT8g02DA88u/tR7A+KrYE2hvSQ9XPWqh5n1lVfbaUwNr2jI5ZAdldMymPm0famtzEiRZ+/KSkMH+qFaD59W6fzP9QwgfghoSCmttOnmgEiakitv4S1QmhkLS66ebaxeYGac8A9bl6E+s3QPZuE0CBBpw5afvScALiOgIeWub9LdVKtYslhz667tNqQB8TZzcSigKZhNbfuu/bHwsPnmiMKapvX8HoHBTAMEj1zXOm2CmTz8+tzAItpY5bQA3EJAQ8qpikMhLRLFtYg3/f00a6tXM7E5xzRe3tHUWZ9njpvXy5w85VrTcPUJUcHsz7e2sT9OOq6W96n7Xv54qZk4s+QQRDi2IwH3EdDghN80OydqsGB54Xrz/sp/em9A+fnlsq3H3m5WLCmuvlA4azv7AtPghzamwt7ifxyGB7NY9RuuGj99nflwUmjIJugMmixZ/ZP3BMBVBDQ4geJaJIK2yv1y2RumP2cW/rTSTGvzvllXO7h+QsMCd/13XloGM9EF6M99uMhs37nbhrMhd7X1vlPSynWF3hMAVxHQ4AyKaxEvGjBRuWzOmOttuaw6zXyVa203fS6qZHr8Mvos2eZtP5u//H2h7UnLX5FeQyoKZ/e9+L29FF2u6NbITnIGDQrkr2QAB3AdAQ3O6FK7te1Qi/RUPqtoKBsFMQUyBbM/znu3xCRwvco17K0WBT1fMgNb9TW3/FeePVemFbPIuzJVYnvzX6abwTGKXl3Q+95vbb3FO+OWm+kLN9tw5g8AnHZCbXvnpuQGXPmkEKetUJ+qM+59YTZn0wCHENDgFIprUR76+0NbmNrK1JamtjZ92j5/tf1tNpjpXtjDDwmtKFWrUtFOZp5z4jHmnitaBt6bmePoVUgKZNrGFN2tefeQWUXhTHdwPnR1K/sseQ2reU8lfThxpf2qgYIxk9eYKfM2cTYNcAgBDU6huBYHYvrmxebyKU/ZclkNAWgYwKfrxnRl2JyznzbXNO6639stIi8418Xp4Vc/uWTqvjAVpH7tw82DVx3rvYX0/GXdwGGB2Uu22BW44V8st++6YeC2i5rbZwCpR0DLAH+Y+445Y9ID5siPLzcVRl5of7QYe6v9fH9eXjq26OdE/kgFimtRFuPXzbKt/yqXHb5iUok6ljOPbmM+OvF+M+3Mv9h7Xct67ZjPD2qv3dfRVKro5j8iVwQc8lc4e/yG46IuQNe7AmdptCX6pxuP994AuICAlsYUrhqM/rV5eN7fzYT1s0ts62jLR59rdaE0326c7z2VpG2hVKG4FrFou/usrx+0P9T6H05h7MtTH7E/giaCD1RZwpmujBo3da3ZvWev90ny+KtiGgJQAAsKZz6tjOmCdIU4n36eVs1uvzi3xJYoADdU2LuP94w0ouCllYOyWNhtiMmtGrxVo5U3hbtIfRt2Me906u+9JZ8OeussUThtfS7u/vwBr4ggvWl1TMFMh/61pRlOfy9oW3xAi4vslmayvTF6mRk2psCuul3ZLSdq2CCTaDtUFAa1wigaKvibN0ihoNfb0S1hIB2xgpaGFF78cHZ67ePNV10etSGs/uE17WeRtB0Uy3ebg6fUgqYpk4niWiiY6dJynS/TdmZ4ONNBf6206nyZ/kUiFeFMq2cfTAwNr6xct91OQv76ye/sgftUrKgByCwEtDSkqoDzjuloHjr2v204U0jTCln9ysEBLZbIbdFw2i5KJYprs5cO+mv1NGf09ebaaU+XmOBVMNMk5uIeL5gX292a0q14rZRdeFp9Ow3qI6gBiBcCWpoaddID5o+tLvfeQn7cdWDt4JM3LvSeStJKXKwt0WSiuDa7FJXL7gtmkeWyuq/1D8deZlb1etX+C4o6zVKtSuVDzVU9G5th93c0/XrkENQAxBUBLYPE6gqrWSm4B0lhJ0iH6m6M2lNcmx38ctlmn99kz5mt2/mj952wctkeoXJZBTXXKJiVJahNKOXycgCIREDLEJroDKJpyIvqn+S9lTR/a6ioMlLHGrneU+pRXJu59N/hLTOGFpXLhteoaOtSW5jaytSWZuRUr4vCg9o1vRqbWkdW8r4TCmrrNu3w3gBg/whoGSJWXUaH6rHD1nebF3lPJZ1d5wTvKfUors08Ouyvs2U6/D90yeioclkd+tfhfw0B+K3/6URB7cruOWbYAx3NLX2a2aCmH66W3gJwEwEtQ3y3KXgas0ON4O3K0lbcNHTgCk1yUlybGfQvEX65rKYzI8tl1fqvcllVvGRClYp61C46vYENauooc7X0FoCb+CdGhoi1GnZa7eO8p5JmbVnmPZVU2opbqlBcm97Grp1hi2VPnjAgqly2V932ReWyqZ4cThQFs9yG+z87V7hjt9m5a4/3BiDbEdAyQHnOny2Icf4s1opbKulguEJapMGLRpVYhYFbFMa0Wtb9m4Eluvi0OqYwptWyT09+yK6ewZj3xq8w/R6dakZMWElQA0BAywRaoQhSnvNnsVbcUo3i2vTgl8vqLtjIcln996dLy3W+TNuZqSiXdZVferthy04zZORighoAAlommFKO82ertm/03kqKteKWahTXus0vl1Uwi1UuW9DzJfNq+9tSWi7rKk15HlW1uJ4jPKi99XmBDXAAsgsBLQPEqpyItRr2wapvvaeSXOk/i4XiWvdoUOPR+e/Zqgx1mS35aY33neJyWQUzV8plXdUyp5p5ZUAHe8el7vX0Kai99tky0++xqfbeT4IakD24LD3NjdgXti7+95PeW0l7+3zgPRXT6tn105/13krSOS91T7ms81d3RxXs6jyTtsyQPCqXHZz/ka3JiJymVRDT5eXaznSxWNZ1unVg3NS15q2xBXZlLZwqPHS9lKZDwwtxk0EXo0tO3SqmXV51+wwgcVhBS3Mbd271nkoKKvbM37baPDj3be8tWqtqDb0nd1Fcm1r6fbat/2NuNE8sGBFVLvv0CTcUlcsSzspHd3z26Fw3cEVNK2jDxhTY66OSrfcp9ewPwhmQHAS0DHVkxSreU4jCWa9/PmzPnumuzSDhf6Dqr+84vr9doXMJxbWpMevHZUXlsjprFlkuq7NlOvz/22bnpmW5rItiBTVKb4HsQEDLUApi/gF6fT1t0n129UOrHEFnucS/jeAPc98x7cb/zk56Xv3dYBvWXEFxbXJpClPTmL/48o6octmTarYsKpfVdmYmlMu6KDKo3di7KaW3QBbgDFoGOPLjy21x6/5oCODvne82K7avN2dMesD7NDZtk77e4Q7nJjsVxHSxdmQgG3T81eb3eRd4bzgY6i3TxeXh/WW+bnXamvtbXkJ/GQAkEAEtA2iF7PezX/PegkUOACigTVg/23uLdt4xHc3gX1xvcqu6uZWic1Daagunrc/F3Z9nJecgqFz2sfnvRw1iiIYxNJVJf5n7Xvhoidm87WdzZbecEmfYAKQPAlqG0HTmGwXjS4QurZh1qJFr+uWcEXi/5g3TnzMf/zClqBNN25+d9v31tzTr5dR9nEFU56DOrcibBHTRtu5yRNnp9/DNgq/MU/kj7VmzcAq7v9r394+mMtNhiAShag71p6nkVtujZ3esQ1AD0hABDWnr0smDzPsrv/HeQjrVyDOTzxjkvaE0Oug/dPFndsAicgpWB/1v3hfUdd6v6RF1vU+RDnQDgUpuwxHUgPRDQEPa0lCDLuCONLHL46ZL7dbeGyLp7J6GKrQ1rj6zcJrk1STmb5qfS7FsGpuRv9kMG11gv4ZTUDuz/dHm0jMalukCdwCpQ0BDWlNA86dPfRTXBlMY04qZLpmPHLBQGLsjt7e5uWlP+ssySKygJqe2qW369cghqAGOIqAhrWmLU1ud4XRuSp1c3PkYokvldb5M4Sy8v0z0e6SL6LWdSX9Z5pq1+Efz/viV5utZ671Pij150/GmQ0tWSwHXENCQ1nTAXcMC4XdAiprsdf9jNtO5Mk1kagAgcphCB/518F8DAEy9Zo/8FdvsTQR+UNN5NPWraesTgFsIaEh7qttQ7UY4dbgV9HgpK7frVC6rDjNVZkRSRYaqMrQNjOzlB7VT2tSyJbgA3ENAQ9qjuDZEpbJaMRu7dob3STGVyiqYUS4LAOmBgIaMkM3FtVop0xmzSevneJ8U00qZWv9VPwKUx7ipa83RNSqZtrlckg4kEwENGSHbimv1n3P48kk2mGlLM5xfLqs7V9sc1dj7FDhwWwt3mX6PTbVfFdD69cwhqAFJQkBDxsiG4lpNYerScnWYBZXL6kqvO3LPZ4IVcfHG6GX2rFo4ghqQHAQ0ZIxMLq7V+ToFsycXjAgsl1Uw65/Xh3JZxNXKddvNW2ML7Dbn7j0l/6ggqAGJRUBDRsm04tp1O380zyz6JGa57E1Ne9pKEcplkUgENSD5CGjIKJlSXKtVMq2WBZXLavhB58sol0WyKai999UKM+bfa+xl7OEU0O7r19LUOrKS9wmAg0FAQ0ZJ9+Jav1x2+PKJga3/msikXBaptmHLTvPuFyvMqG9WFwU1Sm+B+CKgIeOkY3GtJjEVzFSZETmJSrksXBUe1O64JJfSWyCOCGjIOOlUXKtyWVVljFo9xfukmEpldR1Tr7rtvU8ANymoVa96GKtnQBwR0JCRXC+u1UqZDv4roEXSSpnOmKX75CkQaX7BVpPbsCpBDigDAhoykovFtfq1KJhpKzOoXLZvoy42mGlLE8g0WmXr9+hUeyvBld1yzNkd6xDUgFIQ0JCxXCmu1WF/HfpXMAsql72mcVfbYUa5LDLZkJGLzYgJK7230FBBn1Prm/NPqWcqVTzE+xSAj4CGjJXq4loFM9VkxCqXVTDTGTPKZZENNEigWwm0khZOtRyXdW1IUAMiENCQ0VJRXKvhBJ1/CyqXVTC7o/n5lMsiK6mSQ0FNk58ENaB0BDRktGQW1/rlsi8t/dxs3bXd+zREq2RaLaNcFihbUOvTpT5n1JDVCGjIaMkortW5sicXfGDeLBgfs1y2b6PTCGZAhFhBjdJbgICGLJCo4lpNYqrDbPjySYHlsgpm2k6l9R8oXWRQu7tvC0pvkfUIaMh48S6unbR+jg1mqsyIpHJZVWWcX6+T9wmAslJQ04XsCmesniHbEdCQFeJRXPvZmmn2jFlQuawCmYKZAhqAxNu9Zy8hDhmNgIasUN7iWv31uobpj/PeDSyX1RamtjIplwWSR9ugtzw1w5x38jHmotMbmGpVKnrfATIHAQ1Z40CKaxXM3iz4Kma5rA79K5hRLgskX3jprcKZtkQ1+akJUCBTENCQNcpSXOuXy+qM2fLC9fYzn4KZajIolwVSR1ubv3vmP2bO0i3eJyHqTlOHGkENmYKAhqwSq7j21fa32TNqzy8ZHdj6T7ks4JavZ623NxPkryg5/OMHtUvObGjqVCeoIX0R0JBVgoprD6lwiKlWsbL58edC75MQv1z2+ibdbS0HAPfECmoaINCF7LqYXb1qQLohoCGrxCquDadzZbq8XHdlUi4LpIfJczfaHrUZ+Zu9T0Ja5lQzz97Z1nsD0gcBDVlFk5g3zRhi/r1xgfdJMU1iqiqjb6MulMsCaUoBbdjogqKgNvDaVubUNrXtM5BOCGjIClM2LbQTmUHlshoQUDDTWTQAmUEBbfz0deaOi3O9T4D0QkBDRlOprDrMgsple9Vtb8+YUS4L4EBRlItEI6AhI2mlLKhcVrRS9odjL6NcFoBZuW67eeGjJaZfjxyT27BsU9r6OW+MWWbuuaKl9wkQfwQ0ZAy/XFbXMc3dusL7NERnyn6VcwblsgBKGDR8gRkzOTQ0pLNqZQlq/s/57YXNTZ8u/PMEiUFAQ9rzy2UHLxoVNZ3pl8vqjJnu3gQAn66MuuLhKXa7MlxpQU2rZ79+8ruiLc6/3X6CnRQF4o2AhrS16edt5pnFn5hnF30Ss1xW4YzWfwCxqD9NPWrqU4sUFNSe+WCRGTmp+Pq3Y2pVNkPvasd9oIg7AhrSjsKYQpnCmUJaOIUxdZipXJbWfwBlVZagVvOow0y/R6eanbv2eN8J6dyqpnn8huO8NyA+CGhIG9q+1DamtjO1rRlO58ruyD3fBjPKZQGUl4Lae1+tMOOnrYva+tRq2Q8bdnhvJV13XhPTt2sj7w04eAQ0OE8H/nXwXwMAGgQI1+aoxvZ8mQYAKJcFEC86a/bW2AIzburaqKAWROfRnrz5eNM2t7r3CXBwCGhwlioyVJURVC7bqUaencikXBZAIh1IUKt1ZCUzpH9b+xU4WAQ0OKe0clmVyqrDjHJZAMmiac+gs2dBOrSsYc+jUWKLg0VAgzO0UqatzG83zvc+KUa5LIBUGTJysRkxYaX3tn9Xds8x1/Rq7L0B5UNAQ0r55bI6/B/Z+u+Xy+o6plbVGnqfAkDybC3cZfo9NtV+PRBaRdN0J1BeBDSkhKYwX1r6uRmcP8os3FbcKSR+uax6zJoeUdf7FACS743Ry2z9xoGqXvUwM6R/O1OnOufRUD4ENCSVessUzJ5aODKwXPbmpj3NHbm9KZcFkHLlXT3z6YYB3TTAeTSUBwENSaEw5l/HFFQu+5vm55rfNjuXclkAzlDn2Yz8zd5bMHWjlaZpvSPsahpwoAhoSCgFMx38DyqX1falfx0T5bIAABQjoCEhdK7ssfnvB5bL6sC/Dv5TLgsAQDACGuLKL5cdtXpKVDBTRYaqMiiXBQCgdAQ0xIVKZbWV+dmaad4nxVQqq9b/bnXaep8AAIDSENBwUFQu+1T+SDNp/Rzvk2JaKdNW5kk1W3qfAACAsiCg4YBp6/L9lf+0K2axymV1+J/WfwAAyoeAhjLTFOabBeP3BbMPAstlr2/S3dyRe77Jq1rf+xQAAJQHAQ37tXXXdlsuqxWzoHJZBbP+eX0olwUAIE4IaIhJhbJ/zf8osFz26EpHmd80O8fcmdubclkAAOKMgIYopZXLapVMB/8plwUAIHEIaCjil8sOXz4xKpjpXJmqMiiXBQAg8QhosJOYWjHTZGasctnz63UimAEAkCQEtCymclmdL1OXWSSVy2ors1fd9t4nAAAgWQhoWUjXMKlcVgEtkspl1WGmgAYAAFKDgJYltHWplTKdMQsql72kwcl2xYxyWQAAUo+AluF02F+H/hXMgsplf5Vz5r5gdiHlsgAAOISAlqEUzFSTEVQuW63i4bZcVitmlMsCAOAeAlqG8ctln138qVm380fv0xAVyup8GeWyAAC4jYCWIfxy2deWfRHV+k+5LAAA6YWAluZ0ruyphSNtMAsql9X5Mp0zI5gBAJA+CGhpSpOYqsoYvnxSYLmsVsw0mUm5LAAA6YeAlma+3TjfbmXGKpfVGTN1mQEAgPRFQEsTY9fOsFUZQeWyuoapf24fymUBAMgQBDTHaaXsj/PeDSyX1UqZLjCnXBYAgMxCQHOQzpS9WfBVYLmsgtmvcs6wwYxyWQAAMhMBzSF+uawO/y8vXO99GqIpTNVkUC4LAEDmI6A5wC+XfX7J6KjWf79c9rfNzzVHVzrK+xQAAGQyAloKKYwN3hfMhu4LZrHKZa9p3JXWfwAAsgwBLQV0rkxXMWk7M6hctn9eHxvMKJcFACA7EdCSTIEsZ/T1UfdkahJTW5kaAKBcFgCA7HaI9xVJolWx3zQ7x3sz5qSaLc0Hv7zHTDvzL3bVjHAGAABYQUsBnTfr991fKZcFAACBCGgAAACOYYsTAADAMQQ0AAAAxxDQAAAAHENAAwAAcAwBDQAAwDEENAAAAMcQ0AAAABxDQAMAAHAMAQ0AAMAxBDQAAADHENAAAAAcQ0ADAABwDAENAADAMQQ0AAAAxxDQAAAAHENAAwAAcAwBDQAAwDEENAAAAMcQ0AAAABxDQAMAAHAMAQ0AAMAxBDQAAADHENAAAAAcQ0ADAABwDAENAADAMQQ0AAAAxxDQAAAAHENAAwAAcAwBDQAAwDEENAAAAMcQ0AAAABxDQAMAAHAMAQ0AAMAxBDQAAADHENAAAAAcQ0ADAABwDAENAADAMQQ0AAAAxxDQAAAAHENAAwAAcAwBDQAAwDEENAAAAMcQ0AAAABxDQAMAAHAMAQ0AAMAxBDQAAADHENAAAAAcQ0ADAABwDAENAADAMQQ0AAAAxxDQAAAAHENAAwAAcAwBDQAAwDEENAAAAMcQ0AAAABxDQAMAAHAMAQ0AAMAxBDQAAADHENAAAAAcQ0ADAABwDAENAADAMQQ0AAAAxxDQAAAAHENAAwAAcAwBDQAAwDEENAAAAMcQ0AAAABxDQAMAAHAMAQ0AAMAxBDQAAADHENAAAAAcQ0ADAABwDAENAADAMQQ0AAAAxxDQAAAAHENAAwAAcAwBDQAAwDEENAAAAMcQ0AAAABxDQAMAAHAMAQ0AAMAxBDQAAADHENAAAAAcQ0ADAABwDAENAADAMQQ0AAAAxxDQAAAAHENAAwAAcIox/x9btF62UAjmAAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image(\"pic1.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrix of the Householder transformation is\n",
    "$$\n",
    "\\mathbf{H} = \\mathbf{1} - 2 \\mathbf{u} \\mathbf{u}^T.\n",
    "$$\n",
    "Note that $H$ is symmetric and orthogonal.\n",
    "\n",
    "Given two equal-length vectors $\\mathbf{x}$ and $\\mathbf{y}$, a rotation which brings $\\mathbf{x}$ to $\\mathbf{y}$ is a Householder transform with\n",
    "\n",
    "$$\n",
    "\\mathbf{u} = \\frac{\\mathbf{x} - \\mathbf{y}}{||\\mathbf{x} - \\mathbf{y}||}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAFlCAYAAAAAkiT3AAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAAFiUAABYlAUlSJPAAADaaSURBVHhe7d0HuBXVuTfwRVWkN+lNDoKK9BMEKWLkiyZixHuJ4AfRoKio3ICAgoIcBCkCQj6NDWvIEyDkUwwm6MUrVYEA0ouUANI7RzjS4fJfe80+s2fWPmWfmdkzs/+/PD7sGSXKgTPvrPWWVejKVYKIiMiksPqRiIgoisGBiIhsGByIiMiGwYGIiGwYHIiIyIbBgYiIbBgciIjIhsGBiIhsGByIiMiGwYGIiGwYHIiIyIbBgYiIbBgciIjIhsGBiIhsGByIiMiGwYGIiGwYHIiIyIbBgYiIbBgciIjIhsGBiIhsGByIiMiGwYGIiGwYHIiIyIbBgYiIbBgciIjIhsGBiIhsGByIiMiGwYGIiGwYHIiIyIbBgYiIbBgciIjIhsGBiIhsGByIiMiGwYGIiGwYHIiIyIbBgYiIbBgciIjIhsGBiIhsGByIiMiGwYGIiGwYHIiIyIbBgYiIbBgciIjIhsGBiIhsGByIiMiGwYGIiGwYHIiIyIbBgYiIbBgciIjIhsGBiIhsGByIiMiGwYGIiGwYHIiIHLT3zDHx+cGV6iq4GBwCbOHeHfIPIhH5w7ITW0X6wkGix6pJYuXJ7epuMDE4BNSyTcfFe/O2iXrznhBdlr8i31QuXrmk/i4Ree293fNE+8UviIPnTorTF8+KbismBPp7stCVq9RnCojMrAvi0fGrxamz58XyO94WmRey5P2q15QTj9S+UzxR9xei7nXXy3tE5C4EgAHrPxBv7PynuhP5Xvy09VBxW/kb1Z3gYXAIoFEfbxGL1kW2k87+n6/FkrPfyc9md1/fXDxa5y5xf7XWomihIuouETnp6Pkf5QphwdEN6o4QrcqliTm3vSgDRJAxOAQMggKCg6Fws63iv8v8XV3ZGasJBIq0ktXUXSIqqA0//iC3dHf9dFjdEfJ77a2mT4prCxdTd4KLwSFAjmSeF49PWC1On7mo7ghx4bpTYmGrd9RVzu6o1FhuOWE1EYY/vETJ8rf934rfrX5d5hYAq/OxN/cUg9Lul9dhwOAQIMPf3ywT0VabOv1J7L+U/faSm0rFy0RXE41K1VB3iSgvMrbMECO/n6muhChXrKSY3mqg3MoNEwaHgJi7/JB47a/60rhDzReI9aUTq6vmaoIob7BKQImquYcBL1fIL4Rxy5bBIQAOHT8n+kxcLc6c05fFHam4Xay9Zba6SgzefnrW7CgDReMytdVdIoLtWQdE13+Nk3kGw71VW4npLQeKUkWvVXfChcEhAAa9uUGs3ZGpruwuFbkgFt7+urh89X9OQPkdgkT3mu25mqCU99WRtbIi6aQqGYdhN3YTo256SF2FE4ODz32yaL9467Od6iq+Ta3/KvZfk/1W4wRjNTEw7dfsm6CUNHH7bDF005+jzWxYJXzYvJ/4z+pt5XWYsUPax/YcPiPe/8dudZWzaw/WVJ+cgzelvWePyQQ2USo5e/mCrEYavPHjaGDAC9LidmNSIjAAg4NPXbp8RUyYsU2cv5i3raKKJ+qpT855su4vxKz0waHdUyXSwfiLTkuGiY9++FrdiRRuLO0wXjQr6/z3mV8xOPjUrAX7xObdp9RV7sqcqiKuuVJCXRXciIYPymYedldTKsHgvObzB8gfDXhJmtc2I/Adz/nF4OBDO/ZliY/m5i9/UOhKYVH2WMGrjBAMEBQyGnVXd4hSA1YKWDFg5QD4Xpja7KmUfUlicPAZuZ00c5v8Mb8qnCxYcMD2EbaR8KZElCrk4LwNH8gcA3INgFXC/NtHicfqdJbXqYjBwWewYsDKIREVj9dVn/IPSee5t70km+GIUgUG592zdJSYsmOOuhMZnLei40TRruJN6k5qYnDwEeQYkGtIVImzZcV1Z8qrq/zB8jnVvxkotaChrc2iIbKPwdCzVke5YqhZoqK6k7oYHHwCVUmoTkpkO8mswonEtpaGbpoW3WslCrvZB5aLNoufl53PEBmc10tMa9Gf1XkKg4NPvDtnl+xrKKhES1q3nN4nuix7JTplkiisMDgPozCMP+sIBpiPNKTBA/KaIhgcfACjMT5bEnmDKagKJ+vIyqXc6MZi4MzbbiteDfTRhkTxIBggKJgnqmJwHvILYZuo6gQGhyTDMD1sJzmlyKViouypnCdE9q/fRWy76y1t3fYXh1fLqg2iMMGBPNhGwnaSAQEBgYFj6/UYHJIMc5MwddVJFU7UUZ/sJtzysJjcuLdMuM1rO1K7v/rnPQtlDoIoDJBwTl84OGaiKgbnYSuJ+YX4GBySCCe7HTpxTlQuW1zdcYYu74CEG5Jt5pOqMJob5au6Bp9x2z6JKe8jCiIc+o9SVZSsAoIBvg8wUTUVG9vyg1NZfQJVSvuPnhX7jp6RK4n9x87K64PHIz/mdcYSXCl0WSxs80dxsWhkRYJviE9/NkTcVbmpvLbCkYcYSWyFb57prZ5NmUFjFB7Im/VZ82bMfCSslue0fjGl5iMVBINDAGAy64yv96qriFrXlxAlrikiA0lmVqSr02zdzXPE4Urfy7wCls9o7MkJVgnoErVCgEDdN3sgKChQkt11+diY+Uj48zsr/bmUm49UEAwOATDq4y1i0bpj6ipiyEM3ip+3rCw/nz5zMbLKOHF1xXF15YHPX55bLJYVXyI+TX9R3HnDDfKfyw3GE2N+vRVWHkzcURCg4g4l2eaeHYzAeKvpE9xGyicGhwB48rU1tpEa/++/moib6pRWV3bGN0d+35R6rJwkZuxboq6y4f8HAYKdo+RX2ELqu/bt6HwkBIPXm/ThrLAEMSEdALpqpioVrlGf9PAwT2QJPa1lf21uAsGm87cjYo5KJPID5Bew6rUOzsOYbQaGxDE4+BzyCdg2MkOuoUJpZyucDHjbQvJal6OQXdTL2UVN/oGXFVQjmbdDkXDGwTw4oIcSx+Dgc8gfWFWv5G5tNnIMSGKnlbQ30y05tlm+obGLmpINLyvoXzAPzuteo508ypNnnhccg4PPHdRsKVUt737jDpblCBC6rSmUvg5Yb69sIvLK5wdXXg0Mg2IG56F3YXqrgWxscwiDg8+h+sjK7ZWDAdVJn7Yeqv1mQ3PR6K2z1BWRd/Dnzry9afTxoOuZnMPg4HPooLaqXtG7N6Pbyt8oprccqC0DHL75LzFNRkRuQjBAsyb+3Bmw9YkqunurtlJ3yCkMDj6nG+OdW6WS0/CNh8OAdNCFiiU+kZswOK/9khfklqYhMjhvAvtvXMLg4HOHjtsT0jUqlVCfvPNI7Tvlnq4VEtM9Vk2SiWoiNyw4ukEmntdk7lR3hJwRhpxYuWIl1R1yGoODj2GeEobzWXmVc7DCnu4z9X6prrJFlvuvyuoRIicht9X524zo4DycQ4LBeZguzI5ndzE4+FgyylhzM/nW3tpBfGiSs44tIEqUMTiv37qp0bJpdOcvbj9GnvNM7mNw8DFMZLVKxpaSGd7W0EWtG8SHssJ7lr7MJjkqELxgdPpmuHhv9zx1J1IYgcRzbgMkyTkMDj6mWzl4nYzWwdIeo491iUDsC+MoRjbJUSIwOA/9C+YcFvJdWDEkMg6GEsfg4GM408HKyzLWnCAROL/daO03LDpWe62aoq6I8gYnEGLFsPdMZAIxVqkYnPdh837MLyQBg4OP6QbuJTvnYIbAgACha5LDZFcMQyPKDVaZOJa213dToluSlYqXEXPbDNcWQJA3GBx8DKfCWXkxOiM/sLWEo0ax1WSFYWg8apRygsF5KGTAsbQGHF+L/oV4JxeSNxgcfMyP1Uo6SE4jSa1b+uN0OXPjEpEBpc9tFj0vvji8Wt0RshJuafvxHJznAwwOPoX+BpwrbVa2ZDE5rtuP8E2N2nOdHitfk41MRAZjcJ65N2ZEwwfFrPTBHJznEwwOPrXn0E/qUzY/rhrM+tfvIoY0eEBdZcOeMgalbfjxB3WHUhkG56GizTw4D9VvGY26y2vyBwYHn9IO3PN5cICxN/fSNinhQYCT5DAjh1IT/gzgGFoMzjNKnTE4D9tIHJznPwwOPqXLN1StEIzlNkoPMRTNCs1NaJIzRiFQ6kB5Kgbnmc8nR8J5aYdxMgFN/sPg4FO6cxxqVU5ud3ReITE9K/052dVqhT1m85YChV9kcN6gmMF52IJEqSpKVsmfGBx86qBmW8kP3dF5JQ9gaT1U20WN7lfUtLOLOvwwAgOD84yZWyh5xspycuPebGzzOQYHnwpKGWtO0CSHAKHrop59YLkcqkbhhMCPwXn4y3gJMJomMQ6D/I/BwYcysy6I02cuqqsIlLBWKF1cXQUHVg6Yu68rT3x715ciY8sMdUVhgVUCVgvWwXmrO03WbjWSPzE4+JBubEaV8sHZUrLCJE3kIHTbCCO/nxnzEKFgQ14B+QVzXwtWCvHmcJF/MTj40D7NllKyR3UXFKqXsNes03ftO3KbiYINg/NQkWQenIfGSPy+68arkL8xOPjQoRPBzzfooP8BfRBW2JNGgppHjQYTfv/Qu2AenIepvahGwnGeFEwMDj6kPcchwNtKZuigRhmjFR4qKHHlUaPBgt83DM5D17OBg/PCgcHBh8JQqZQTbDXojhpFcxya5HjUaDAgkCO/YB6ch05ndDyj85mCjcHBh3QNcEHPOZhhL3p6q2e1R41ivAbGbLBJzt8QEDBR1To4DzOSODgvHBgcfOb8xctyIqtVkBrg8gIBAudA6JrkMKAPg/rYJOdPOHsBW0k4iwEQDDBNlYPzwoXBwWfinf5WpHAhdRUeeKigxLFmiYrqTjaUQmLUN/kHVnNIOuPUNiNw49wFbCPptgkp2BgcfEZ7+ltABu4lArXv89qOlNUtVjgkCIcFUfKhPBXnO6Nc1XBHpcYy8czBeeHE4OAzYU9G68gu6jh71ThmlEeNJhdKjJF4Xnlyu7oj5NnO89pmcHBeiDE4+Iz2HIeK4U/wITmNZildFzVWD+ZRz+QddK9jxWBUkOH3Z2qzp8TrTfpof68oPBgcfEbb4xDibSUz7FtPvrW3uorVa9UU8dWRteqK3IacQt+1b9sG5y1uP0Y8VqezvKZwY3DwGV3OoUYKrBwM2K4YdmM3dZUNDyg0yZnPBCB3oN8Eg/MwGNGA+VgrOk7k4LwUwuDgM9qheyErY83NqJse0o51RrUMmuS2Zx1Qd8hpkcF5g2MG52HsCVYMuqoyCi8GBx9BfwP6HMzKliwmSpUoqq5SB/a1decKY+8bNfbsonYeqsMwOM845xs5BczCmtaiPwfnpSAGBx/RdUaHvVIpHjyYprccqO2iRldutxWvsovaQRic123FhOjXFKXFOIdjSIMH5DWlHgYHH9Elo6um2JaSWaTz9rm4R432WDWJXdQFhGCAbnTz4Dx8vZd2GC/HrFPqYnDwEX2PQ3hmKiUCFTJ4g9UdFPP5wZXyLAhKDHI36F/A19GArTwknnUBmVILg4OPaM9xSKFKpXgw4XNum5e0TXKow8eWCOUPBuch8WwenIcqsU9/NkT7dabUw+DgI7oT4FI152DVrGw9+eDSNV5hS8Rcdkk5m7h9tm1w3vRWA2WVGBvbyMDg4CPxhu5RBA6Pmdayv7qK1W/dVFltQ/GdvXxBDs4bvPHjmMF5i9uNEd1rtJPXRAYGB584feaiyMy6oK4iihctLCqULq6uCPAQw2FBVnjY/W716zxqNA6U/rZf/IJtcB4Sz1iVEVkxOPhEKg7cSxTOJY531Cgqb3jUaKxlJ7aK5vMHxAzOwwgMDM7TJfqJgMHBJw5qBu6F6fQ3p01u3Ft7hgD20XGSHJvkIj764Wu5YrAOzsNfzC9QThgcfELXAJdqYzPyC0eNYmvESp49sGRYSjfJYZsNeRhstRn5hcjZGRkcnEd5wuDgE9ptJZax5ghvvjgHQnfYDLaW7ln2cko2yWFw3j1LR4k3dv5T3YlUe6F/QRdMiXQYHHxCe44Dcw65QhkmTpJD1Y2V7KJe+VpKBQicv43+BfN4czk4rx0H51H+MDj4xJ5DP6lP2cJ8PKiTsF2CJjndqWQob0XpZirAr7XN4udjBuehdwGD89jYRvnF4OADly5fkRNZrbhyyDuMe4jX3YtjRtH4FWYZW2bEDM7D1wFfD93ZGER5weDgA/rT364RRQoXUleUF5jgirdkXRUOVg/mGv+wMMp3R34/U92JBErkF3Qjz4nyisHBB7Snv7GMNSH3V2stzzfWQeUOZgqFBQbnYRvJPDgPk1TR2MbBeVRQDA4+oD39rTzLWBP1ZN1fiBENH1RX2ZCYxjkQ5mawoELCuc2iITIBbRjS4AE5wRZnMRAVFIODD+w/xu5op2U06q6t55fbMMuC3UWNHApKVVGyCsgvYDsNp7axsY2cwuDgAzzHwR1vNX1CbjNZoVu46/KxgeuixuA8bI0N2PBBtDwX5anzbx8ly1WJnMTg4AMHj9uDQ1VuKxUY3qLxRh3vqFEEiKB0USOQoesb4zAM+HUh8dyqXJq6Q+QcBgcf4NA99xglnboELQbSIQfh9yY5Y3AefjRgywwrBg7OI7cwOCQZ+hvOX7ysriLKliwmSpUoqq6ooNAchyY53YMU1Ut91ryprvwHKwWsGMyD81CNxcF55DYGhyQ7pNlS4sA952G8BsZs6Jrk8AAeummauvIHrGaQW0COAbkGQJDD4Lxn6v1SXhO5icEhyfYc1vU4cEvJDRjQh0F9ujfucds+iRlUl0wYO45qJFQlGSKD8yZwcB55hsEhyfQD91ip5BY8XDHqW2fA+g+SftSobnAezq3A4DzdcEEitzA4JJn2HAdWKrkKD1scFmSFrZxeq6Yk7ahRdDqj4xmdzwYMzpuVPli7HUbkJgaHJDuo6Y6udT1XDm7DMaO6o0axv49zILxuksPgPMxIMg/OwxYYB+dRsjA4JJl+6B7fEr2A1UP3Gu3UVTY8oM0VQm7CvwvTVM2D89JKVhNL24/n4DxKKgaHJDp95qLIzIpUohiKFy0sKpctrq7IbdNa9hd3VW6qrrIZTWdIDrsF5y5gG8mc58B/CxLPutPtiLxU6MpV6nOorNmeKWZ+vVfs2J8lTpyKPIDLly4m0huVF4O7N5DXOu//Y7f4csWh6M+5oXpJ8c7AZvKz03bsyxJPvrZGXUXUrXqdmDq4uboiL+Dtvf2SF8SazJ3qTjZ0IaME9trCxdQdZyDhjFPqjPlIMCjtfjH25p7sXyBfCOXKAQ/4wW9tECu/Pxl9yAM+//eKw+Llj7eoO7GGvrtRzLgaUMw/58w597pn92kG7vH0N+9hfx9NctjOsUJyGklqJ7uoUTJrHpyHwPNh835iwi0PMzCQb4QuOCAw4AGfk8XrjsmVhRkCA4KJVVoN98Yf6yqVODYjOdA9jXHXui5qbPugzLWgEGDQ1NZv3dRosMG/b3H7MeKR2nfKayK/CFVwmPPtwWhguO/2amLepNtFq4b62TN//ya7XBABRRcYwM1Dd7Q9DhUZHJIF85dmpT+nLRvF2z4a5RKFHEb7xS/EDM67rfyNYnWnyRycR74UquCAElAEg1vqlhb9HrhB3nu6a+RHq+37IonGBWuORgPKtcWLiP/6j/piQt/sLtSWcYKLE+IdD0rJgxzD9JYDtds7GLFhfrjnFQ4XSl84KGZwHlYKWDFwcB75VWgT0mZPTFoj/r3fXnXy4ZAW4tk/rpc5BgSGgQ+miTuaVVJ/1309X1lpOwXug+dbsM/BB97bPU87kA9BA1Ne81pminOr+6z5Y3Q+En4+cgu6HgsiP0mJUtZ4eYMXpm6KJp+9DgyXLl/RHg/KnIM/YCQ2upOtkCvosWpSzCpAB//c4I0fi17fTYkZnDe3zXAGBgqElAgOjWqXVp9iHVDVQt3vrOlpYADdlhL6G4oULqSuKNnQnYzzqK1Q+oqDguJ1UaM3AkeRTtw+W92JDP1b2mGctqeCyI9SIjh0aVtVfbJDjuLRX9VRV97Rnf7GgXv+g7MTMIvJCglmBABrFzUCBgbn4ZwIA44qRcezrlSWyK9SIjhANU0VEJrixj5+i7ryFk9/CwbkCNCDoDtqFAPyECCMeUgYnIfEs3lw3oiGD8ocBQfnUdCkTHCocDUQWJ05F3sCm5f2axrgGBz8yRiCpztqFJVIXf81Try8daZtcB6CQkaj7vKaKGhSJziUsc8rOnv+kth7xN6I5gVdMrpKeQYHvypXrKQco6ErPcUojBGbZ6ir7MF52E4iCqqUCA7T/2ev7IrWWb0ttlPaK/s03dE8Ac7fapaoKEZcXQkUKhS/aAAJZySeOTiPgi70wQFNbn/5Kv44jS0/nFKfvMWcQ/CgQxqjL3JqDbqnSgtZskoUdKEODtgy+uCfu+X2EZrcdPYlYVvp+Knz4vzF2HxHqRJF5V/kP8ZMJHRI4zPyCfESzM9v/FPSjxolckKog8PEGduivQyPd6krx29b6RLDbuOqITjkuQ7fDI+OzUDOYf7to2SyWTdiQzbJrXwtaUeNEjkltMEBU1Y37opsGbVvUlH2OugqltAhbZ3Q6jZdcGC+wX9w2H+bRc9HH/TmQXnILaAHQgcBAkeN4ucTBVUogkP/19eJzgO/kec0YCvJPH4b/Q0vPdxIfk6rUUr+aGVOVmOyK2YxYVKrW3SrFVYq+Ytx2D9Oa4OetTqK+e1Gx1QroXta10ENKGnt/O0IsfeMvhCCyO9CERyMhy0e8r8b9100MCDPMKbPzfIzxJuwunjdURlUkLx+d84uOaRv/b/dW03ou6MZHPwCiWf0LuABj60jzFia1qK/9jQ4rB7ijcTAlhQChPm0N6KgCHxwwJaQ+eQ2gzFltWbl7JEUzdLKavMO+PkIKq9M+14mr7HaGJTDUaIFxZyDP+kSz8gtYMZSPAges9IHaxvkAOM0jEBDFCSBDw6bd9tLUXMav93j5zXVJz0ED6w2zEHFafrgwLlKyWQ9jKfuddfLRra8jOZGg9ynrYfKH3WQs0DQQcAhCorAB4efzl6SM5IAQQGD9EY9elPcKau4jwN9cCCQAT8P15jO+s7AZq4GBpxJnZkVu9IpXrSwnMhKydN37dvRMdxYBazoOCFfjWyRU+QGayuYAOWt6JEgCoqUOOzHT3bsyxJPvrZGXUXgcB8c8kPJg8Rxk/m/FycuZMkHPM6Tvvv65urv5t3bu76UgSYe5C9y2qYi8otQJKSDZJ+mUsnNc6opbz7a87X4UeUFIr0Kk+Ke15CTnCqYYPjmv8hT5oj8jisHj81asE9WRJk90KG66PvreuqKvIREcd91b8vjPK2wVbS0w/i4uYR4ZJ/D0lFyIJ+OkcROxmC+93d/Jf60Z4H4LnNHNEmOQYEIaAPTfi2vrXZkHZQVXDP2LZY/B4n6e6u0EtNbDVT/RPiN2DJdljd/l/lvdUeIFmVvkDmpkY16qDt21b/sLQ6cPaGusuHnrrpjkrryJ64cPKbvcbhGfSIvYSup/ZIXYgLDHZUai2fr3yc/Y+XQbcWEfCeSc6tgwv8fjg/1uov63mWjxWNr/igWHdsYUz2F8ycGbfxITNr+mbqTDYHh7qUvy9WO8XPMPzfsPjmwTDT46inx8vd/jQkMgGvcx9c1nsXtxogOFe1nxuDn4vfBzxgcPMYyVn9YcHSDPJhnTeZOdSeyJTSvbYaY1Ph30bd6vP0nkkjG23Wj0vEr4/CA7bbi1YS2rhKBB9g/Dq1SV3qTdsQGByMwmA8vMqBLPOzw8H74uz9of/1m+LpiRaZTv2RV8dtad6irbNWuLa8NGn7C4OAxXQNc1QoMDl7CW3DnbzOiR3ziTX9qs6fEW02fjFYboemtWdnIVh+SzPnJExjzmGYfWC6v0VX988pN5Gcz/HP3XH34Gv8dbumz5k35AEPAeq/Z02L7XW9pjyzF9od59fD79e/FfTCmlw9/cHh09R9lEMfXamG70eL//+z5uAMXsVWXH7+qknuJdLIxOHjo0uUrXDkkEbZz8KDEX8ZWkTFI77E6neW1AQ8BnP5mjMvou/YdudrIjXUeE96wMY/p89uGaY8axXgOBAg3t2owEwp73Pg1PlrnLvk2Gy9pjpPtYMCGD6IrDePh+FLD38jrILz1OuGhmu3lr/395k/LX+8D1W4THeP8unPaItpwyj5jq1etjuqTfzEh7SEEhofHxi7t0d/wl5fS1RW5BW/nqEAyP+Dx4EYHNA7xiQe9D52WDBNnL1+Q5zTgIB/dWzcgYdlj1aTogx7zmKZefVM3xm7gPraydFtJyHVgS8tYuXih0Gdd1adsCCIvNuwm/uNf4+U1fq1ftHlJBhQScmWF/IwOVmUIvlbWLb1fVWkpXxb8jisHD+lOf6vCLSXXIa+At3lzYJCD9K6uGHIKDIC3bjzgATOSuizLPifabPTWWTFnSOvmMWE1Yh3eZ8B/W69VU9SVNxAIrJAoxZkUgBUCA0OseBVdoFshgDWR3fXqCiQIGBw8pNtSQgMcuWfGviWyIsmYroo387E395IP7nj7x1YIJEbjWqSC6dXothRWFKg8Qv8CGNtR8Rrd5DbW1QChK4/Ffyu2c7zSolx99SkW8gz4dbzR5HEGBg1dUIWDmpJVVDuZS1nxc3WrCz/itpKH0N+APgezR+6uLf5v51rqipyCh/fILTPlG70BD7xZ6c8l1PkMGKBnJJmHNHhA/L5+F9F1+djo2A3MY0JgyMvYDeQkMLEVwcVqcuPeov/V/2+3ocIGpa06yC/kVL+fqJz+nU7Dg9iNXgLkrHQFCrp/n3VLKd7Wkx9x5eAhDtzzBrZ28CA3B4bIvKSJCQcGMFcwoSns1q9/Hw0MyBmgYS6v85iQnJ7Wsr82x4DVgxdHjeK/WQd74m4EhrDAVqPO1qz96lMESoEXmhLV+LoGJTAAg4OHdDkHngDnLGz7IOmL5LABAQGBIV5TWl4ZW0ZGzsA4pwFVQEgm63IJOfnP6m3F5Ft7q6tYOGo0L9VRBYEtI93W2ilNToWyxQuq1lwUXiDM955rYC8A8DNuK3novheWyamsZn97+WeibEn7ITKUf18cXi0rkk5eyFJ3hNz7H9HoQUergLBawHhvmHDLwwXeAsL5EXiQWOHBjQ5bY7XiBnT/WnsZUKG07a431RXp6L5uYGwbYdXQbMGAaHDACwR6aYKEKwePHD913hYYSpUoysDgEDxcUUlkBAY8WLENhKohp8tDsa3wYfN+cnKrE7kBJMgfqX2nusqGBwt6IIxkutMQSHUPuHiNb5StYanq6lMs488fGgiNwICqL+SogobBwSOHjp9Tn7Kx+a3gjGoh4/Q2QHkq3rhRZeQW/H8XJH9hhbdK3cFCbnVRY5AcqqPiiTcOgiIaxAkO2NbE186chB5Y/9eBrPpicPCILhldlQP3CkQOzltsH5yH/IKbWzFuwOpmesuB2mSnUT5r3dNOFMorX9vxd/lZl3OAeDX7FNG+YvbZ9GZbT+8Xw7dEypoBndU59Ub4GYODR3TJaFYqJQ77/kg8G+MeINHEsF/gQY3jRnWJc5S+YoVkrI4Shb3wZ9a9KwMN/n1r7pis/k6sbVcfchQfRmnoYIyG0deAr++om4Jb9cXg4BHdqG5uKyUGNeZYMRhbLXjrxtA8bM04nV/wGgIbchm6AIceC8x4SpQxZdV4eGU07C63O7AnbmXt6iW7eM1wBrysBHkGFYODR3TbSjzHIX/w1ozx2dbBeVgt5HT6WtCgWggBQrflg8CYsWWGuspZ6X/0kPOTMA8Ib7Tm8dvda7SLbnfoHnIIINh+MuBrjoYu5iKy3Rgn7wD4mqKZMcgYHDyiS0jXqnKd+kS5QU8Bxmy/sfOf6o6QeQU0nsWrOw8yYyigbiU08vuZcox4ThAMjBwFBsV1XDIsGhjw4DKf4tYyzhiNP+z4XP6I5DWCEpKszEVkq6pZcQGC+l/TB6ur4GJw8ABKWFHKalakcCE5kZVyh8F56QsHxzSF4c0XFUkYWRFWd1VuKktmdbCCyqmL+n+OrFOfYmFVYn1w/bZWJ+0qBQEGKw+cdgbo8A3627CTGpfWd8Mb23VBx+DgAf3YDOYb8gIPQOvgPPQu4M1X90ALG5TMotHOCttqv1v9etyjRlE1YxVv/DausT+eEwTjIIyZ9pJuFAYCaFCrk6wYHDxw8IQ9ONRgpVKuMOkUZzgb2yMIBthqiTfxNKwGpd2vbbbD1wUzpHTnQ4ARPJFwxsM/p/HbWBFg2B4CiAE/Dw87dP2at6EowpyTAXzt/nDrY+oq+Dg+wwOYxIqJrGa/bldNPNM152qHVIWHHg7NMc9HMpK0BZ2PFFRYKWDekm4rCU1/6O0IaglvEKHyCytac9nqxy1+H7fENYi4cvCAtoy1IreVdJA0bbP4ec3gvAkpGxhANsm1elZ71CiaATH+21hhkft6r34j5pyGZ+vfF6rAAAwOHtBVKlWpwDJWq6+OrJWJZ5zDbMCWClYMusNxUg0CxNzbXtKOBcfX7J5lLxe4SY5yh5Je85nRYR1xzuDgAf2obuYczCZuny3uWToqOrgMx2ticB6SsUFvbHMSti/mtR2pPd4UyWlsPVHBVf+yt6zUwnBCM5T1mucmoSw4rIl6BgcPsFopPgzOQ9XN4I0fR9965eC89u4OzguySOPfSFGpeBl1JxtyEl4eNRpGyCcYW0Y4rAfXWCmgV8Qo6wVdWXCYMDi4TBcYKpQuLooX5Zce4y8wBuOjH75WdyLjsJFcRRMYxYf8Cyq3dOW8U3bMkSsxSoy5nwZBIu2rvjIwmLeS4pUFhwmfUC7TD9zjqgGD85rPH2AbnIcVA6tu8gbJaTTJ6bbdsBIzT6ulvDOOfo0nFQIDMDi4jOc42GGlYB2c93qTPqEYnOc1HDWKr50OtuuQ5Kf80TUQGjBILxUCAzA4uIxlrNmQU8DYBzy0jPwC9s3nthkunqn3S3lN+Yehg7rGQHyN0SRnXp1R7lqUu0GuDgzmZsCF7UanRGAANsG5LOPDLeKbDcfUVcQLPRuKTs0rqavUgMF5qKQxv8miJBMH9od5PpKXMDkVA/KssE2H7TrzA48oN1w5uEw3OqNqivU4oAa/zaIhMYEB2yFL249nYHDQW02fEPdXa62usmH7DudrO33UKIUbg4PLUr2MFQfUoOPZfGj9iIYPilnpg7WVNpQ45GvQG6Lrosb8pa7Lx7KLmvKMwcFFmVkX5Lhus1IlioqyJYupq3DDoTTY8zYPzsM2Ukaj7vKanIev8az057SjRlCFg5lV7KKmvGBwcJH29LcU2FIypoXiUBoD9ruxjXRv1VbqDrkFOYa5bV7SlgRjZhVyE0S5YXBwkS441Ah5pZIxOA/bSQYcWrO0wzjtTCByB3I5CBC6rTuUEmMcOlFOGBxcdPC4Lt8Q3plKSDgj8WwenIdzCFCqqhv1QO7CMarYxtP1jozeOivmyFUiKwYHF+05Yu+OrlI+nNtKGNmAwXkoWQUMzkP3Lg6RYWNb8uB87Wkt+6urWAPWf5DjUaOU2hgcXJQK3dFIbqKpDcPejEQn9rrntxstHql9p7ym5MIRn7qzn/H71WvVlLhHjVJqY3BwkTbnEKJtpXiD81Z3mix/JP/A9p7uqFFMxe2y/JW4R41S6mJwcAlKWI+fOq+uIooULhSaaiWMZMDgPPOQMqwUsGLg4Dx/wuoBqwgrnKHRackwNslRDAYHlxw6Ed4tJd3gPDx4kGNAroH8C/kH5CGs8HuJAMEmOTIwOLhEN6q7aoVgBwfsUWMUNHIM2I4AHN+JaiTdlgX5DwI5KphQyWSFrSUcNWr83lJqY3BwSdjGZmDrAdVI5kNk0LeAg//Rx0DBgd4H9EDo5lohOY0kNbuoicHBJbptpaCWsaJvAQf/mwfnodMZHc+c9BlMOXVRo7wVK0RKbQwOLgnLygHjFnSD87A1wcF5wYb5S5jDpPt9RN/KuG2fqCtKRQwOLtHlHIJWxorBeShzNA/OwzRVDs4LD0xwxSRXXaPi0E3TYsqUKbUwOLhE1wAXlG0lBINuKybEDM7D/jS2kXAOA4ULzoDAWRA6GNL3xeHV6opSCYODC7CldOly7AF7FUoXFyWu8f8YiV0/HRbtl7wQM1YBpY9IPHNwXng9Vqez3C60QmK624pXcz10n8KHwcEF+oF7/t+fX3B0g0w8r8ncqe4IebbzvLYZHJyXArBdiPOoreQI9uVj2UWdYhgcXKBLRvu9MxoTOjt/mxEdnIc96KnNnhKvN+nDwXkpBL/fuq1DNMkhQLCLOnUwOLhAV8Zaq7I/k9HYNsC+cr91U6O17ShvxIH02Gqg1IIXAXS6xztqFGdRs4s6NTA4uGB/QLqj8RaIMRjv7Z6n7gjRqlyaWNFxIgfnpTBUpX36syHao0YxUws5CDbJhR+Dgwv2HfN/zgHf5OkLB8UkGnvW6ihXDDVLVFR3KFUhxzSv7UhtkxyqlzBChcKNwcEFfj/H4c97FsoVw94zx+Q1thIm3PKwrHfn4Dwy4CUBAULXJIc/Q+iDoPBicHBYZtYFcfrMRXUVgRLWsiWT/9DFVgC+oXt9NyVmcN6c214Ug9Lul9dEZihfnnvbS9qiBHRQo5OawonBwWF+HZuBwXlIJppHImBPeWmH8eLu65urO0R2SE5Pb/WsuoqFGUw8ajScGBwcdlCzpVSjYnKDA6pM2ix6PqbTFYPzkHjWJR2JrFDeGu+o0R4rX+NRoyHE4OAwXaVSlSRWKmFwHhLP5gamYTd2k9UoHJxH+YEzO3TbjwgQOAeCTXLhwuDgsP26SqUkrRxGb50luv5rXMzgvOmtBopRNz3ExjZKCAoXUNVmhT9jOEnOKHKg4GNwcJgfcg74Ru2xcpIYvvkv0Xp0DM5b3G6M9gxhovxAk5zugCf0zXT+doTMb1HwMTg4TBccvGyAMwbnzdi3RN2JDM5D4ll3NCRRfmHViW1JNExayS5q05h3Ci4GBwedv3hZHD91Xl1l82rlgMF5SDybB+dhBAYG5+mamYgShS1KlEDrTgJEchpNcuyiDjYGBwclc0vp7V1fysF5xmA0vN1hcB7+Yn6B3IAXDgSIeEeNDlj/gbqiIGJwcFAyTn/D2xkG5/Vd+3b0TQ3frFgtcHAeuQ2l0J+2HqqtfMOkXxRFUDAxODhIe/qbi6O6IwnAjJjBecgroH8BeQYiL2BI4/SWA7UrVBRFmP98UnAwODjIyzJW5BXQv4A8g0EOzmvHwXnkPTRVYgtTp+/ad2S/DQULg4ODvMo5yMF5S2IH56F3AYPz2NhGyfJI7Tvln0MrbHf2WDWJXdQBw+DgIN3xoFXLO/ewxjcZlukYnGdubENSEF3PRMmGP4c4WtYKf15xDgS7qIODwcFBbq4c8M2FwXnmBB+SgcgvcHAe+cnkW3vHPWoUf4aNijryNwYHhyAZfenyFXUVUaF0cTmuu6DwtoX8gnlwHgICGts4OI/8Btuc01r21x41uj3rgLhn6ctskgsABgeH6MpYnahUQkCwDs4b0uABuZWEsxiI/AiHRs1p/aL25QXFFJj5xSY5f2NwcMihE86f/oYtJPOB7sgvIOk89uZebGwj38PLy/x2o7VNcl8dWSt6rZqirsiPGBwcos03JFjGimCApLN5cB7KU+ffPko7EZPIrxAYECB0VXSY/4XDgsifGBwcojvHIZHuaJSndvpmuCxXNWDvFoln3aAzIr/D1hKOGtWdTz5x+2weNepTDA4OOejAthLqwJFfWHlyu7oTGZyHFYNuaU4UFHjBQZJatx06YMMHPGrUhxgcHKLbVspPQhojBrBiMA/Oe71JHw7Oo9BAeSsOC9LBUaPmbn9KPgYHB2RmXRCnz1xUVxEoYUUpa26QU8DQPAzPM/ILlYqXkYPzdM1EREGGo0ZRbWeFP/s4B2LDjz+oO5RsDA4O0A3cy8uWkjE4D+O2DZHBeRM4OI9CC9V2GLVhhUIMnCSHA6so+RgcHLBPs6WU29gM1HrjYB7zUhrLbgzOw5GeRGGyvV8/cXrNGnUl5HaprrMfL0xokjt6/kd1h5KFwcEBukqlnFYOKOHD4DzzGxIGls1KH8zBeRRKxz7/XKxq3lxs7NpVBgnk0WalPyfHfVuh4RNNcuyiTi4GBwdoG+Di9DigdwGH/5sb29BJysF5lAqOzp4dDRKXlq2UBwXpuqhRuYdeH3ZRJw+DgwP2HM59dAaCARJu5sF5OH8X/QuYhU+UShAk1rRvLw7d3VXMPH+HtlR79oHlot+6qeqKvMbg4IBDmlHd5gY4DBtD/4L5wJO7KjeViWcOzqNUdnLBAnHil93ErFc2izab7XkGFGtkbJmhrshLDA4FdP7iZXEk87y6ymasHCKD8wbHDM4blHa/mNtmOAfnESmXliwXY4d+Lf4w7BvRbMNRdTdi5PczYyr6yBsMDgUU7wyHIoULydEAGJx38kKWvI/xAR827ycbgdjYRmTXdMMRMWXYEvmXOUhgewnbTOSdQleuUp8pAcs2HRfD3489/rBJw1LicKv5MfORsKeKMducj0SpaHm9euLsrl3qKu/WNK4k/talvljSupos3sCMJt05EeQ8BocC+mTRfvHWZzvVlRDnip8Su1r/U+wptEfdEbJcD1UZnI9EqSrR4GDYXq+s+Kh7I7GlfUOxuP0Y5uo8wG2lAtp/LHtbKbPMAbG8xbSYwIBOUPxhZmAgSlzazkwxeuxyMe7pz8SLk3qzSc4DXDkUELaUsLW0r+o6saXBV+JKocvyPnIKyC1gloyTLp0+LS4cjU3YucXLf9fls2fF+YMH1ZW7rly8KM7t3auu3Hdu9271yX0FeTvPL3wN8bXMi/z8s3lRsllTUXdEhqh0//3qDjmNwaGAeo//TlSf+6rovHCmukNEXinVrJmoM2IEg4QLuK1UQM/+Jk10TKuprojISxjFgW7rTd26OboyIQaHAmtcr4xoULOUuiIiL5VISxMNP/xQ3DR9uihUtKi6S05gcCCiwDGCQvrmzaLqI48wMLiAOQcH7MrIELtHjlRX7ipSqpQoVqmSunJX4WuvFcWrVlVX7vLy34UHSfEa3pVCXlu3rvrkPi//XdfUrJnnh/LmHj3Ej8uWqavEISjUfvFFUaVnTwYElzE4EJHr1nbqJOcoJYpBwXvcViIi38Jq8sapU7l9lAQMDkTkOwgK9SdPFq137hTVHnuMQSEJGByIyDfMQaFm//4yF0XJweBAREnHoOA/DA5ElDQMCv7F4EBEnitarpy4YcIEBgUfY3AgIs8gKGAWEoJCrUGDGBR8jH0OROQ6zD8q2bSpXCUgQJD/MTgQkeswFI/lqMHC4EBERDbMORARkQ2DAxER2TA4EBGRDYMDERHZMDgQEZENgwMREdkwOBARkQ2DAxER2TA4EBGRDYMDERHZMDgQEZENgwMREdkwOJCv7X/zTbGmY0exsFCh6F+4xv2cLK9bN+bnmH8uEeWOU1nJlzIXLRLbnnlGZK1fr+7YVenZUzSaNk1dxTq/f7/Y+sQT4tjnn6s72VqtXStKNmmirohIhysH8p2sdevEhvvuyzEwwKE//1kcmzNHXcUqXr26qPb44+oq27V16jAwEOUBgwP5zuarK4KLmZmi5K23yrf8ZgsXiqJly6q/G+vAu++qT3lz/UMPqU9ElBMGB/Kd6k8+KQPDTVdXBnjLL9uhg6jYpYv6u7F020aGrI0b1ads13fvrj4RUU4YHMh3qj/1lGi1bl3M9k/5X/xCfbKLt7X0kyU4IEfBLSWivGFwoEDAgz0e3QoBMhcvVp8iKv/mN+oTEeWGwYECA9tLOuf27lWfsqHa6ezu3eoq8nPjbU0RkR2DAwVG2dtvV59i6aqaDkydqj5F1Bo0SH0iorxgcKDAKNOmjfoUK2vtWvUpAj0O5jwEtqS4aiDKHwYHCozSLVuqT7FQ9mq27403Yu7VGjxYfSKivGJwoMBAYxtKXHWMlQJWDebRGrWHDmWFElECGBwoUEo1bao+xTJWCv9+/vnoZ3RD13jmGfmZiPKHwYEC5bpbblGfYv20aZNcPWCkhqHu6NFytUFE+cfgQIFStm1b9SkWKpa29+unroSoeO+9OfZGEFHOOJWVAgejt3OCOUzNFi1iroGoALhyoMCJ1wxnqDVkCAMDUQExOFDgxKtYAgSO2leDAxEVDIMDBc41NWuqT7GwnXTz9OnqiogKgsGBAqdknIqltDfeYHUSkUOYkKZAsialczoylIjyjysHChxMXDVDDuKG8ePVFRE5gcGBAgXjMbb89rfqKpJnaMDtJCLHMThQoGx94omYcxrqjRmTa2krEeUfgwMFxpZevWLOjEaeAUeKEpHzGBzIV5bXrSuTzdssA/MwadU8NwmrBSagidzD4EC+gXyCsWV0/OoKAddZ69aJDV26iG1PPy3vAxLQ7GcgchdLWck3MFV1w333qSs9BIYmX3zBBDSRy7hyIN/4celS9UmPgYHIOwwO5BsYux0PRnAzMBB5h8GBfKPM7bfHDNXDSW6oSGr897+LxnPmMDAQeYg5ByIisuHKgYiIbBgciIjIhsGBiIhsGByIiMiGwYGIiGwYHIiIyIbBgYiIbBgciIjIhsGBiIhsGByIiMiGwYGIiGwYHIiIyIbBgYiIbBgciIjIhsGBiIhsGByIiMiGwYGIiGwYHIiIyIbBgYiIbBgciIjIhsGBiIhsGByIiMiGwYGIiGwYHIiIyIbBgYiIbBgciIjIhsGBiIhsGByIiMiGwYGIiGwYHIiIyIbBgYiIbBgciIjIQoj/BTZb4NZCJSBYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Image(\"pic2.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write a function which rotates a given vector, \n",
    "$$\\mathbf{x} = \n",
    "      \\begin{pmatrix}\n",
    "            x_1 \\\\\n",
    "            x_2 \\\\\n",
    "            \\vdots \\\\\n",
    "            x_n\n",
    "       \\end{pmatrix} \\;\\; \\mbox{into} \\;\\; \\mathbf{y} = \n",
    "      \\begin{pmatrix}\n",
    "            ||\\mathbf{x}|| \\\\\n",
    "            0 \\\\\n",
    "            \\vdots \\\\\n",
    "            0\n",
    "       \\end{pmatrix}$$ \n",
    "using a Householder transformation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def householder(vec, sign = 1):\n",
    "    \"\"\"Construct a Householder reflection to zero out 2nd and further components of a vector.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    vec : array-like of floats, shape (n,)\n",
    "        Input vector\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    outvec : array of floats, shape (n,)\n",
    "        Transformed vector, with ``outvec[1:]==0`` and ``|outvec| == |vec|``\n",
    "    H : array of floats, shape (n, n)\n",
    "        Orthogonal matrix of the Householder reflection\n",
    "    \"\"\"\n",
    "    #vec = np.asarray(vec, dtype=float)\n",
    "    vec = vec.copy()\n",
    "    if vec.ndim != 1:\n",
    "        raise ValueError(\"vec.ndim = %s, expected 1\" % vec.ndim)\n",
    "    outvec = vec.copy()#np.asarray(vec.copy(), dtype=float)\n",
    "    outvec[0] = sign*np.sqrt(np.sum(np.square(vec)))\n",
    "    outvec[1:] = 0\n",
    "    u = np.true_divide((vec - outvec),(np.sqrt(np.sum(np.square(vec - outvec)))))    \n",
    "    H = np.eye(vec.shape[0], dtype = float) - 2*u*u.reshape(-1,1)\n",
    "    return outvec, H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test your function using tests below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "ename": "AssertionError",
     "evalue": "\nNot equal to tolerance rtol=1e-07, atol=0\n\n(mismatch 33.33333333333333%)\n x: array([ 3.741657e+00,  0.000000e+00, -1.110223e-16])\n y: array([3.741657, 0.      , 0.      ])",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-5-d000d15424d2>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mv1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mh\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mhouseholder\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msign\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0massert_allclose\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mh\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mv1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mv\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 6\u001b[1;33m \u001b[0massert_allclose\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mh\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mv\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mv1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32mC:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Anaconda3_64\\lib\\site-packages\\numpy\\testing\\nose_tools\\utils.py\u001b[0m in \u001b[0;36massert_allclose\u001b[1;34m(actual, desired, rtol, atol, equal_nan, err_msg, verbose)\u001b[0m\n\u001b[0;32m   1394\u001b[0m     \u001b[0mheader\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m'Not equal to tolerance rtol=%g, atol=%g'\u001b[0m \u001b[1;33m%\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mrtol\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0matol\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1395\u001b[0m     assert_array_compare(compare, actual, desired, err_msg=str(err_msg),\n\u001b[1;32m-> 1396\u001b[1;33m                          verbose=verbose, header=header, equal_nan=equal_nan)\n\u001b[0m\u001b[0;32m   1397\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1398\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Anaconda3_64\\lib\\site-packages\\numpy\\testing\\nose_tools\\utils.py\u001b[0m in \u001b[0;36massert_array_compare\u001b[1;34m(comparison, x, y, err_msg, verbose, header, precision, equal_nan, equal_inf)\u001b[0m\n\u001b[0;32m    777\u001b[0m                                 \u001b[0mverbose\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mverbose\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mheader\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mheader\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    778\u001b[0m                                 names=('x', 'y'), precision=precision)\n\u001b[1;32m--> 779\u001b[1;33m             \u001b[1;32mraise\u001b[0m \u001b[0mAssertionError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    780\u001b[0m     \u001b[1;32mexcept\u001b[0m \u001b[0mValueError\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    781\u001b[0m         \u001b[1;32mimport\u001b[0m \u001b[0mtraceback\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mAssertionError\u001b[0m: \nNot equal to tolerance rtol=1e-07, atol=0\n\n(mismatch 33.33333333333333%)\n x: array([ 3.741657e+00,  0.000000e+00, -1.110223e-16])\n y: array([3.741657, 0.      , 0.      ])"
     ]
    }
   ],
   "source": [
    "# Test I.1\n",
    "\n",
    "v = np.array([1., 2., 3.])\n",
    "v1, h = householder(v, np.sign(v[0]))\n",
    "assert_allclose(np.dot(h, v1), v)\n",
    "assert_allclose(np.dot(h, v), v1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test I.2\n",
    "\n",
    "rndm = np.random.RandomState(1234)\n",
    "\n",
    "vec = rndm.uniform(size=7)\n",
    "v1, h = householder(vec, np.sign(vec[0]))\n",
    "\n",
    "assert_allclose(np.dot(h, v1), vec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.000000000000001"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(h[i][i] for i in range(len(h)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See how many tests (there are 2 of them - I.1 and I.2) your function passed. Sum diagonal elements in matrix `h`. You will need this for Google Form."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part II. Compute the $\\mathrm{QR}$ decomposition of a matrix.\n",
    "\n",
    "Given a rectangular $m\\times n$ matrix $\\mathbf{A}$, construct a Householder reflector matrix $\\mathbf{H}_1$ which transforms the first column of $\\mathbf{A}$ (and call the result $\\mathbf{A}^{(1)}$)\n",
    "\n",
    "$$\n",
    "\\mathbf{H}_1 \\mathbf{A} =%\n",
    "\\begin{pmatrix}\n",
    "\\times & \\times & \\times & \\dots & \\times \\\\\n",
    "0      & \\times & \\times & \\dots & \\times \\\\\n",
    "0      & \\times & \\times & \\dots & \\times \\\\\n",
    "&& \\dots&& \\\\\n",
    "0      & \\times & \\times & \\dots & \\times \\\\\n",
    "\\end{pmatrix}%\n",
    "\\equiv \\mathbf{A}^{(1)}\\;.\n",
    "$$\n",
    "\n",
    "Now consider the lower-right submatrix of $\\mathbf{A}^{(1)}$, and construct a Householder reflector which annihilates the second column of $\\mathbf{A}$:\n",
    "\n",
    "$$\n",
    "\\mathbf{H}_2 \\mathbf{A}^{(1)} =%\n",
    "\\begin{pmatrix}\n",
    "\\times & \\times & \\times & \\dots & \\times \\\\\n",
    "0      & \\times & \\times & \\dots & \\times \\\\\n",
    "0      & 0      & \\times & \\dots & \\times \\\\\n",
    "&& \\dots&& \\\\\n",
    "0      & 0      & \\times & \\dots & \\times \\\\\n",
    "\\end{pmatrix}%\n",
    "\\equiv \\mathbf{A}^{(2)} \\;.\n",
    "$$\n",
    "\n",
    "Repeating the process $n-1$ times, we obtain\n",
    "\n",
    "$$\n",
    "\\mathbf{H}_{n-1} \\cdots \\mathbf{H}_2 \\mathbf{H}_1 \\mathbf{A} = \\mathbf{R} \\;,\n",
    "$$\n",
    "\n",
    "with $\\mathbf{R}$ an upper triangular matrix. Since each $\\mathbf{H}_k$ is orthogonal, so is their product. The inverse of an orthogonal matrix is orthogonal. Hence the process generates the $\\mathrm{QR}$ decomposition of $\\mathbf{A}$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write a function, which receives a recangular matrix, $A$, and returns the Q and R factors of the $QR$ factorization of $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qr_decomp(a):\n",
    "    \"\"\"Compute the QR decomposition of a matrix.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    a : ndarray, shape(m, n)\n",
    "        The input matrix\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    q : ndarray, shape(m, m)\n",
    "        The orthogonal matrix\n",
    "    r : ndarray, shape(m, n)\n",
    "        The upper triangular matrix\n",
    "        \n",
    "    Examples\n",
    "    --------\n",
    "    >>> a = np.random.random(size=(3, 5))\n",
    "    >>> q, r = qr_decomp(a)\n",
    "    >>> np.assert_allclose(np.dot(q, r), a)\n",
    "    \n",
    "    \"\"\"\n",
    "    a1 = np.array(a, copy=True, dtype=float)\n",
    "    m, n = a1.shape\n",
    "    r = a1\n",
    "    q = np.eye(m,m)\n",
    "    H = np.zeros((m,m))\n",
    "    for i in range(n):\n",
    "        q_i = np.eye(m,m)\n",
    "        a_i = r[i:,i:].copy()\n",
    "        x = np.array(a1[i:,i])\n",
    "        x_out, h = householder(x, np.sign(a[1][1]))\n",
    "        a1[i:, i:] = h@a_i\n",
    "        q_i[i:m,i:m] = h.copy()\n",
    "        q = q@q_i\n",
    "    print(r)\n",
    "    return q, a1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Might want to turn this on for prettier printing: zeros instead of `1e-16` etc\n",
    "\n",
    "np.set_printoptions(suppress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 7.28010989  6.73066763  6.59330707]\n",
      " [ 0.          0.83553169  3.13888933]\n",
      " [ 0.         -0.          6.90475747]\n",
      " [ 0.         -0.         -0.        ]\n",
      " [ 0.          0.         -0.        ]]\n"
     ]
    }
   ],
   "source": [
    "# Test II.1\n",
    "\n",
    "rndm = np.random.RandomState(1234)\n",
    "a = rndm.uniform(size=(5, 3))\n",
    "a = np.asarray([[4,4,2], [4,4,4], [2,2,8], [1,1,4], [4,3,1]])\n",
    "q, r = qr_decomp(a)\n",
    "\n",
    "# test that Q is indeed orthogonal\n",
    "assert_allclose(np.dot(q, q.T), np.eye(5), atol=1e-10)\n",
    "\n",
    "# test the decomposition itself\n",
    "assert_allclose(np.dot(q, r), a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now compare your decompositions to the library function (which actually wraps the corresponding LAPACK functions)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "ename": "AssertionError",
     "evalue": "\nNot equal to tolerance rtol=1e-07, atol=0\n\n(mismatch 100.0%)\n x: array([[-0.549442, -0.361311,  0.399255,  0.025288, -0.638376],\n       [-0.549442, -0.361311,  0.109599, -0.029503,  0.744772],\n       [-0.274721, -0.180656, -0.814166, -0.44012 , -0.187935],...\n y: array([[ 0.549442,  0.361311, -0.399255,  0.122609,  0.627001],\n       [ 0.549442,  0.361311, -0.109599, -0.143044, -0.731501],\n       [ 0.274721,  0.180656,  0.814166, -0.406205,  0.253026],...",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-21-6103ecb12979>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mqq\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrr\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mqr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0massert_allclose\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mqq\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrr\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0ma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 6\u001b[1;33m \u001b[0massert_allclose\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mqq\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mq\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      7\u001b[0m \u001b[0massert_allclose\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mrr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mr\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Anaconda3_64\\lib\\site-packages\\numpy\\testing\\nose_tools\\utils.py\u001b[0m in \u001b[0;36massert_allclose\u001b[1;34m(actual, desired, rtol, atol, equal_nan, err_msg, verbose)\u001b[0m\n\u001b[0;32m   1394\u001b[0m     \u001b[0mheader\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m'Not equal to tolerance rtol=%g, atol=%g'\u001b[0m \u001b[1;33m%\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mrtol\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0matol\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1395\u001b[0m     assert_array_compare(compare, actual, desired, err_msg=str(err_msg),\n\u001b[1;32m-> 1396\u001b[1;33m                          verbose=verbose, header=header, equal_nan=equal_nan)\n\u001b[0m\u001b[0;32m   1397\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1398\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Anaconda3_64\\lib\\site-packages\\numpy\\testing\\nose_tools\\utils.py\u001b[0m in \u001b[0;36massert_array_compare\u001b[1;34m(comparison, x, y, err_msg, verbose, header, precision, equal_nan, equal_inf)\u001b[0m\n\u001b[0;32m    777\u001b[0m                                 \u001b[0mverbose\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mverbose\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mheader\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mheader\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    778\u001b[0m                                 names=('x', 'y'), precision=precision)\n\u001b[1;32m--> 779\u001b[1;33m             \u001b[1;32mraise\u001b[0m \u001b[0mAssertionError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    780\u001b[0m     \u001b[1;32mexcept\u001b[0m \u001b[0mValueError\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    781\u001b[0m         \u001b[1;32mimport\u001b[0m \u001b[0mtraceback\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mAssertionError\u001b[0m: \nNot equal to tolerance rtol=1e-07, atol=0\n\n(mismatch 100.0%)\n x: array([[-0.549442, -0.361311,  0.399255,  0.025288, -0.638376],\n       [-0.549442, -0.361311,  0.109599, -0.029503,  0.744772],\n       [-0.274721, -0.180656, -0.814166, -0.44012 , -0.187935],...\n y: array([[ 0.549442,  0.361311, -0.399255,  0.122609,  0.627001],\n       [ 0.549442,  0.361311, -0.109599, -0.143044, -0.731501],\n       [ 0.274721,  0.180656,  0.814166, -0.406205,  0.253026],..."
     ]
    }
   ],
   "source": [
    "# Test II.2\n",
    "\n",
    "from scipy.linalg import qr\n",
    "qq, rr = qr(a)\n",
    "assert_allclose(np.dot(qq, rr), a)\n",
    "assert_allclose(qq, q)\n",
    "assert_allclose(rr, r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Did `q` and `r` agree with `qq` and `rr`? If not, explain why. You will need this for Google Form."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part III. Avoid forming Householder matrices explicitly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the special structure of the Householder matrices: A reflector $\\mathbf{H}$ is completely specified by a reflection vector $\\mathbf{u}$. Also note that the computational cost of applying a reflector to a matrix strongly depends on the order of operations:\n",
    "\n",
    "$$\n",
    "\\left( \\mathbf{u} \\mathbf{u}^T \\right) \\mathbf{A}  \\qquad \\textrm{is } O(m^2 n)\\;,\n",
    "$$\n",
    "while\n",
    "$$\n",
    "\\mathbf{u} \\left( \\mathbf{u}^T \\mathbf{A} \\right) \\qquad \\textrm{is } O(mn)\n",
    "$$\n",
    "\n",
    "Thus, it seems to make sense to *avoid* forming the $\\mathbf{H}$ matrices. Instead, one stores the reflection vectors, $\\mathbf{u}$, themselves, and provides a way of multiplying an arbitrary matrix by $\\mathbf{Q}^T$, e.g., as a standalone function (or a class)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write a function which constructs the `QR` decomposition of a matrix *without ever forming the* $\\mathbf{H}$ matrices, and returns the $\\mathbf{R}$ matrix and reflection vectors. \n",
    "\n",
    "Write a second function, which uses reflection vectors to multiply a matrix with $\\mathbf{Q}^T$. Make sure to include enough comments for a marker to follow your implementation, and add tests. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qr_decomp_new(a):\n",
    "    \n",
    "    a1 = np.array(a, copy=True, dtype=float)\n",
    "    m, n = a1.shape\n",
    "    r = a1\n",
    "    q = np.eye(m,m)\n",
    "    H = np.zeros((m,m))\n",
    "    for i in range(n):\n",
    "        q_i = np.eye(m,m)\n",
    "        a_i = r[i:,i:].copy()\n",
    "        x = np.array(a1[i:,i])\n",
    "        outvec = x.copy()#np.asarray(vec.copy(), dtype=float)\n",
    "        outvec[0] = -np.sign(a1[0][0])*np.sqrt(np.sum(np.square(x)))\n",
    "        outvec[1:] = 0\n",
    "        u = np.true_divide((x - outvec),(np.sqrt(np.sum(np.square(x - outvec)))))\n",
    "        uai = u.reshape(1,-1)@a_i\n",
    "        r[i:, i:] = np.eye(len(x))@a_i - 2*(u.reshape(-1,1))@uai\n",
    "        q_i[i:m,i:m] = (np.eye(len(x)) - 2*(u.reshape(-1,1)@(u.reshape(1,-1))))\n",
    "        q = q@q_i   \n",
    "    \n",
    "    \n",
    "    return q, r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test III.1\n",
    "\n",
    "rndm = np.random.RandomState(1234)\n",
    "a = rndm.uniform(size=(5, 3))\n",
    "a = np.asarray([[4,4,2], [4,4,4], [2,2,8], [1,1,4], [4,3,1]])\n",
    "\n",
    "q_new, r_new = qr_decomp_new(a)\n",
    "\n",
    "# test that Q is indeed orthogonal\n",
    "assert_allclose(np.dot(q_new, q_new.T), np.eye(5), atol=1e-10)\n",
    "\n",
    "# test the decomposition itself\n",
    "assert_allclose(np.dot(q_new, r_new), a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And again, compare your decompositions to the library function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "ename": "AssertionError",
     "evalue": "\nNot equal to tolerance rtol=1e-07, atol=0\n\n(mismatch 80.0%)\n x: array([[-0.549442, -0.361311,  0.399255,  0.025288, -0.638376],\n       [-0.549442, -0.361311,  0.109599, -0.029503,  0.744772],\n       [-0.274721, -0.180656, -0.814166, -0.44012 , -0.187935],...\n y: array([[-0.549442,  0.361311, -0.399255, -0.262916,  0.58227 ],\n       [-0.549442,  0.361311, -0.109599,  0.306736, -0.679315],\n       [-0.274721,  0.180656,  0.814166,  0.337478,  0.339313],...",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-24-44daba825260>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;31m# Test III.2\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[0massert_allclose\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mqq\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mq_new\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      4\u001b[0m \u001b[0massert_allclose\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mrr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mr_new\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Anaconda3_64\\lib\\site-packages\\numpy\\testing\\nose_tools\\utils.py\u001b[0m in \u001b[0;36massert_allclose\u001b[1;34m(actual, desired, rtol, atol, equal_nan, err_msg, verbose)\u001b[0m\n\u001b[0;32m   1394\u001b[0m     \u001b[0mheader\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m'Not equal to tolerance rtol=%g, atol=%g'\u001b[0m \u001b[1;33m%\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mrtol\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0matol\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1395\u001b[0m     assert_array_compare(compare, actual, desired, err_msg=str(err_msg),\n\u001b[1;32m-> 1396\u001b[1;33m                          verbose=verbose, header=header, equal_nan=equal_nan)\n\u001b[0m\u001b[0;32m   1397\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1398\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Anaconda3_64\\lib\\site-packages\\numpy\\testing\\nose_tools\\utils.py\u001b[0m in \u001b[0;36massert_array_compare\u001b[1;34m(comparison, x, y, err_msg, verbose, header, precision, equal_nan, equal_inf)\u001b[0m\n\u001b[0;32m    777\u001b[0m                                 \u001b[0mverbose\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mverbose\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mheader\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mheader\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    778\u001b[0m                                 names=('x', 'y'), precision=precision)\n\u001b[1;32m--> 779\u001b[1;33m             \u001b[1;32mraise\u001b[0m \u001b[0mAssertionError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    780\u001b[0m     \u001b[1;32mexcept\u001b[0m \u001b[0mValueError\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    781\u001b[0m         \u001b[1;32mimport\u001b[0m \u001b[0mtraceback\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mAssertionError\u001b[0m: \nNot equal to tolerance rtol=1e-07, atol=0\n\n(mismatch 80.0%)\n x: array([[-0.549442, -0.361311,  0.399255,  0.025288, -0.638376],\n       [-0.549442, -0.361311,  0.109599, -0.029503,  0.744772],\n       [-0.274721, -0.180656, -0.814166, -0.44012 , -0.187935],...\n y: array([[-0.549442,  0.361311, -0.399255, -0.262916,  0.58227 ],\n       [-0.549442,  0.361311, -0.109599,  0.306736, -0.679315],\n       [-0.274721,  0.180656,  0.814166,  0.337478,  0.339313],..."
     ]
    }
   ],
   "source": [
    "# Test III.2\n",
    "\n",
    "assert_allclose(qq, q_new)\n",
    "assert_allclose(rr, r_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sum all elements in matrix `q_new` and `r_new` separately. You will need this for Google Form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.7070537225440219"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q_sum = sum(q_new)\n",
    "sum(q_sum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-9.724906106319466"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r_sum = sum(r_new)\n",
    "sum(r_sum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}