{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gram–Schmidt orthogonalization\n",
    "\n",
    "Chapter 4.4 illustrates a hand technique for computing orthonormal vectors q₁,q₂,… from arbitrary vectors a,b,… with the property that the first k vectors in the original set span the same subspace as the orthonormal set, and this is true for k=1,2,3,...\n",
    "\n",
    "We will move this hand technique to the computer in this notebook.  Some of you will notice that on the computer one can combine operations in a simpler block fashion.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6×4 Matrix{Int64}:\n",
       "  5  10  6  8\n",
       " 10   6  5  1\n",
       "  2   1  7  3\n",
       "  7   5  4  8\n",
       "  1   3  8  3\n",
       "  1   8  9  4"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# start with four arbitrary independent vectors in ℝᵐ\n",
    "# with random entries from 1 to 10.\n",
    "m = 6\n",
    "a₁ = rand(1:10,m)\n",
    "a₂ = rand(1:10,m)\n",
    "a₃ = rand(1:10,m)\n",
    "a₄ = rand(1:10,m)\n",
    "A = [a₁ a₂ a₃ a₄] # show them as the columns of a 6×4 matrix A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6×4 Matrix{Float64}:\n",
       "  5.0   5.61111   -3.16215     1.37162\n",
       " 10.0  -2.77778   -0.0979465  -3.70534\n",
       "  2.0  -0.755556   6.16936     1.22464\n",
       "  7.0  -1.14444   -0.324354    4.20193\n",
       "  1.0   2.12222    5.22283     0.0756904\n",
       "  1.0   7.12222    1.49913    -1.74319"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The vₖ are vectors, but they are all orthogonal and\n",
    "#span([v₁]) = span([a₁])\n",
    "#span([v₁ v₂]) = span([a₁ a₂])\n",
    "#span([v₁ v₂ v₃]) = span([a₁ a₂ a₃] )\n",
    "#span([v₁ v₂ v₃ v₄]) = span([a₁ a₂ a₃ a₄])\n",
    "v₁ = a₁\n",
    "v₂ = a₂ - v₁*(v₁'a₂)/(v₁'v₁)\n",
    "v₃ = a₃ - v₁*(v₁'a₃)/(v₁'v₁) - v₂*(v₂'a₃)/(v₂'v₂)\n",
    "v₄ = a₄ - v₁*(v₁'a₄)/(v₁'v₁) - v₂*(v₂'a₄)/(v₂'v₂) - v₃*(v₃'a₄)/(v₃'v₃)\n",
    "\n",
    "# gather into a matrix V with orthogonal but *not* orthonormal columns\n",
    "V = [v₁ v₂ v₃ v₄]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6×4 Matrix{Float64}:\n",
       " 0.372678    0.571756   -0.358733    0.223061\n",
       " 0.745356   -0.283047   -0.0111116  -0.602583\n",
       " 0.149071   -0.0769889   0.699888    0.199157\n",
       " 0.521749   -0.116616   -0.0367966   0.683342\n",
       " 0.0745356   0.216248    0.592508    0.0123092\n",
       " 0.0745356   0.725734    0.170071   -0.283488"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# now we normalize\n",
    "q₁ = normalize(v₁)\n",
    "q₂ = normalize(v₂)\n",
    "q₃ = normalize(v₃)\n",
    "q₄ = normalize(v₄);\n",
    "\n",
    "# Gather into a matrix Q with orthonormal columns\n",
    "Q = [q₁ q₂ q₃ q₄]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Matrix{Float64}:\n",
       "  1.0          -9.21828e-17   8.99478e-17   1.35877e-16\n",
       " -9.21828e-17   1.0           1.0142e-17    1.42552e-16\n",
       "  8.99478e-17   1.0142e-17    1.0          -2.09229e-16\n",
       "  1.35877e-16   1.42552e-16  -2.09229e-16   1.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#check that Q has orthonormal columns\n",
    "Q'Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Matrix{Float64}:\n",
       "  0.0          -9.21828e-17   8.99478e-17   1.35877e-16\n",
       " -9.21828e-17   0.0           1.0142e-17    1.42552e-16\n",
       "  8.99478e-17   1.0142e-17    0.0          -2.09229e-16\n",
       "  1.35877e-16   1.42552e-16  -2.09229e-16   2.22045e-16"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q'Q - I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q'Q ≈ I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Matrix{Float64}:\n",
       " 180.0          -1.42109e-14   1.33227e-14   1.28786e-14\n",
       "  -1.42109e-14  96.3111        9.32625e-15   6.94211e-15\n",
       "   1.33227e-14   9.32625e-15  77.7003       -1.16249e-14\n",
       "   1.28786e-14   6.94211e-15  -1.16249e-14  37.8113"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compare to what happens if we didn't normalize:\n",
    "V'V # = diagonal matrix (orthogonal columns, but not orthonormal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Matrix{Float64}:\n",
       " 13.4164  11.7766   10.3605   8.86974\n",
       " -0.0      9.81382   9.2715   6.67879\n",
       "  0.0      0.0       8.81478  1.38213\n",
       "  0.0      0.0       0.0      6.14909"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# What does this triangular structure say?\n",
    "round.(Q'A, digits=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QR factorization\n",
    "\n",
    "How do we do all this at once on a computer? We ask the computer to factor the matrix as $QR$ (orthonormal columns times upper triangular)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}}\n",
       "Q factor:\n",
       "6×6 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}}:\n",
       " -0.372678    0.571756    0.358733   -0.223061   -0.0840586  -0.590504\n",
       " -0.745356   -0.283047    0.0111116   0.602583    0.0189054  -0.0272178\n",
       " -0.149071   -0.0769889  -0.699888   -0.199157   -0.619405   -0.242243\n",
       " -0.521749   -0.116616    0.0367966  -0.683342    0.131189    0.478182\n",
       " -0.0745356   0.216248   -0.592508   -0.0123092   0.74464    -0.204877\n",
       " -0.0745356   0.725734   -0.170071    0.283488   -0.192913    0.566789\n",
       "R factor:\n",
       "4×4 Matrix{Float64}:\n",
       " -13.4164  -11.7766   -10.3604   -8.86974\n",
       "   0.0       9.81382    9.2715    6.67879\n",
       "   0.0       0.0       -8.81478  -1.38213\n",
       "   0.0       0.0        0.0      -6.14909"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F = qr(A)  # returns a \"factorization object\" that stores both Q (implicitly) and R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Matrix{Float64}:\n",
       " -13.4164  -11.7766   -10.3604   -8.86974\n",
       "   0.0       9.81382    9.2715    6.67879\n",
       "   0.0       0.0       -8.81478  -1.38213\n",
       "   0.0       0.0        0.0      -6.14909"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R = F.R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6×4 Matrix{Float64}:\n",
       " -0.372678    0.571756    0.358733   -0.223061\n",
       " -0.745356   -0.283047    0.0111116   0.602583\n",
       " -0.149071   -0.0769889  -0.699888   -0.199157\n",
       " -0.521749   -0.116616    0.0367966  -0.683342\n",
       " -0.0745356   0.216248   -0.592508   -0.0123092\n",
       " -0.0745356   0.725734   -0.170071    0.283488"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q2 = Matrix(F.Q)  # extract the \"thin\" QR factor you would get from Gram–Schmidt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Matrix{Float64}:\n",
       " -1.0  -0.0  -0.0   0.0\n",
       "  0.0   1.0   0.0  -0.0\n",
       " -0.0   0.0  -1.0   0.0\n",
       " -0.0  -0.0  -0.0  -1.0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "round.(Q'Q2, digits=5) # almost I, up to signs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Matrix{Float64}:\n",
       " -13.4164  -11.7766   -10.3604   -8.86974\n",
       "   0.0       9.81382    9.2715    6.67879\n",
       "   0.0       0.0       -8.81478  -1.38213\n",
       "   0.0       0.0        0.0      -6.14909"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R # Recognize this matrix?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q2*R ≈ A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element Vector{Float64}:\n",
       " 0.5450097629781777\n",
       " 0.8599801580391012\n",
       " 0.7036387925825908\n",
       " 0.7553540639899048\n",
       " 0.7234080262946185\n",
       " 0.14725528162868073"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = rand(6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Vector{Float64}:\n",
       "  0.08995466028300682\n",
       " -0.0768023071730735\n",
       "  0.07513430689992019\n",
       "  0.03688677949256471"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A \\ b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Vector{Float64}:\n",
       "  0.08995466028300675\n",
       " -0.07680230717307351\n",
       "  0.07513430689992023\n",
       "  0.03688677949256475"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv(A'A) * A'b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Vector{Float64}:\n",
       "  0.0899546602830068\n",
       " -0.07680230717307351\n",
       "  0.07513430689992022\n",
       "  0.03688677949256465"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R \\ (Q2'b)[1:4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Vector{Float64}:\n",
       "  0.08995466028300678\n",
       " -0.07680230717307349\n",
       "  0.07513430689992023\n",
       "  0.03688677949256459"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F \\ b  # the factorization object F can be used directly for a least-square solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.7.1",
   "language": "julia",
   "name": "julia-1.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}