QR_Factorization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QR Factorization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Chapter 1 we saw that the LU factorization essentially captured the elimination process and stored the result in a way that allowed us to use elimination to solve similar systems without having to carry out the elimination again.  The QR factorization accomplishes something similar for the orthogonalization process.  Given a matrix $A$ with linearly independent columns, the QR factorization of $A$ is a pair of matrices $Q$ and $R$ such that $Q$ is orthogonal, $R$ is upper triangular, and $QR=A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Structure of orthogonalization \n",
    "\n",
    "If $A$ is an $m\\times n$ matrix with linearly independent columns, it must be that $m \\ge n$.  The matrix $Q$ then will be $m\\times n$ with orthonormal columns, and $R$ will be $n\\times n$ and upper triangular.  For example, if $A$ is a $6\\times 4$ matrix, the matrices have the following structures, with the $A_i$ and $U_i$ being vectors in $\\mathbb{R}^6$.\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "A =  \\left[ \\begin{array}{c|c|c|c} & & & \\\\\n",
    "A_1 & A_2 & A_3 & A_4 \\\\ & & & \\end{array} \\right] \\hspace{2cm}\n",
    "Q =  \\left[ \\begin{array}{c|c|c|c} & & & \\\\\n",
    "U_1 & U_2 & U_3 & U_4 \\\\ & & & \\end{array} \\right] \\hspace{2cm}\n",
    "R = \\left[ \\begin{array}{cccc} * & * & * & * \\\\ 0 & * & * & * \\\\ 0 & 0 & * & * \\\\ 0 & 0 & 0 & *  \\end{array}\\right]\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "The columns of $Q$ are the result of applying the orthogonalization process to the columns of $A$.  If we suppose that this is the case, let's explain why $R$ must be triangular by looking at the product $QR$ one column at a time.  For the first column we have the following vector equation which specifies the linear combination of the $U$ vectors that form $A_1$.\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\left[ \\begin{array}{c|c|c|c} & & & \\\\\n",
    "U_1 & U_2 & U_3 & U_4 \\\\ & & & \\end{array} \\right]\n",
    "\\left[ \\begin{array}{c} r_{11} \\\\ r_{21} \\\\ r_{31} \\\\ r_{41} \\end{array} \\right]\n",
    "= r_{11}U_1 + r_{21}U_2 + r_{31}U_3 + r_{41}U_4 = A_1\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know however that $U_1$ is the unit vector in the direction of $A_1$.  This means that $r_{21}=r_{31}=r_{41}=0$ and \n",
    "$r_{11} = |A_1|$.  Let's also note that  $|A_1| = U_1\\cdot A_1$.\n",
    "\n",
    "For the second column we have a similar equation.\n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\left[ \\begin{array}{c|c|c|c} & & & \\\\\n",
    "U_1 & U_2 & U_3 & U_4 \\\\ & & & \\end{array} \\right]\n",
    "\\left[ \\begin{array}{c} r_{12} \\\\ r_{22} \\\\ r_{32} \\\\ r_{42} \\end{array} \\right]\n",
    "= r_{12}U_1 + r_{22}U_2 + r_{32}U_3 + r_{42}U_4 = A_2\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "We know from the orthogonalization process that $U_2$ is built by subtracting from $A_2$ the component that is in the $U_1$ direction.  Thus, $A_2$ is a linear combination of $U_1$ and $U_2$.  This means that $r_{32}=r_{42}=0$ and $r_{12}$ and $r_{22}$ are the coordinates of $A_2$ with respect to $U_1$ and $U_2$, which we can compute as $r_{12} = U_1\\cdot A_2$ and \n",
    "$r_{22} = U_2\\cdot A_2$.\n",
    "\n",
    "Carrying out the same reasoning for the last two columns, we find that in general $r_{ij} = U_i\\cdot A_j$ and that $r_{ij} = 0$ for $i>j$ because the span of $\\{U_1, U_2, ..., U_i\\}$ is equal to the span of $\\{A_1, A_2, ..., A_i\\}$.\n",
    "\n",
    "To arrive at this conclusion this more succinctly, we could multiply the equation $A=QR$ by $Q^T$, which gives $Q^TA=Q^TQR$ and $Q^TA = R$ since $Q^TQ = I$.  If we then understand the matrix product $Q^TA$ as a collection of dot products between rows of $Q^T$ and columns of $A$, we have again that $r_{ij} = U_i \\cdot A_j$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 1:  Finding a QR factorization\n",
    "\n",
    "Let's compute the QR factorization for a specific $6\\times 4 $ matrix.\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "A = \\left[ \\begin{array}{rrrr} \n",
    "1 & 3 & 2 & 0 \\\\ \n",
    "0 & -1 & -1 & 0 \\\\ \n",
    "2 & 2 & 1 & 1 \\\\\n",
    "-1 & 1 & 3 & 4 \\\\\n",
    "-4 & 0 & 1 & -2 \\\\\n",
    "0 & -1 & -2 & 5 \\\\\n",
    "\\end{array}\\right]\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Q\n",
      "[[ 0.21320072  0.71960864 -0.32643805  0.03625532]\n",
      " [ 0.         -0.2638565   0.01525411 -0.00954087]\n",
      " [ 0.42640143  0.38379127 -0.10982963  0.09922508]\n",
      " [-0.21320072  0.33581737  0.74745161  0.52474801]\n",
      " [-0.85280287  0.28784346 -0.32338723 -0.07251063]\n",
      " [ 0.         -0.2638565  -0.46677591  0.84150499]] \n",
      "\n",
      "QTQ\n",
      "[[ 1.  0. -0.  0.]\n",
      " [ 0.  1. -0.  0.]\n",
      " [-0. -0.  1.  0.]\n",
      " [ 0.  0.  0.  1.]] \n",
      "\n",
      "Q\n",
      "[[ 0.21320072  0.71960864 -0.32643805  0.03625532]\n",
      " [ 0.         -0.2638565   0.01525411 -0.00954087]\n",
      " [ 0.42640143  0.38379127 -0.10982963  0.09922508]\n",
      " [-0.21320072  0.33581737  0.74745161  0.52474801]\n",
      " [-0.85280287  0.28784346 -0.32338723 -0.07251063]\n",
      " [ 0.         -0.2638565  -0.46677591  0.84150499]] \n",
      "\n",
      "R\n",
      "[[ 4.69041576  1.2792043  -0.63960215  1.2792043 ]\n",
      " [ 0.          3.78993883  3.90987361 -0.16790868]\n",
      " [-0.         -0.          2.07455958  1.19287176]\n",
      " [ 0.          0.          0.          6.55076331]] \n",
      "\n",
      "QR\n",
      "[[ 1.  3.  2. -0.]\n",
      " [-0. -1. -1. -0.]\n",
      " [ 2.  2.  1.  1.]\n",
      " [-1.  1.  3.  4.]\n",
      " [-4.  0.  1. -2.]\n",
      " [ 0. -1. -2.  5.]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from laguide import DotProduct as Dot\n",
    "from laguide import Magnitude\n",
    "\n",
    "# Define A\n",
    "A = np.array([[1, 3, 2, 0], [0, -1, -1, 0], [2, 2, 1, 1], [-1, 1, 3, 4],[-4, 0, 1, -2],[0, -1, -2, 5]])\n",
    "\n",
    "# Slice out the columns of A for processing\n",
    "A_1 = A[:,0:1]\n",
    "A_2 = A[:,1:2]\n",
    "A_3 = A[:,2:3]\n",
    "A_4 = A[:,3:4]\n",
    "\n",
    "# Carry out Gram-Schmidt process\n",
    "U_1 = A_1/Magnitude(A_1)\n",
    "W_2 = A_2 - Dot(A_2,U_1)*U_1\n",
    "U_2 = W_2/Magnitude(W_2)\n",
    "W_3 = A_3 - Dot(A_3,U_1)*U_1 - Dot(A_3,U_2)*U_2\n",
    "U_3 = W_3/Magnitude(W_3)\n",
    "W_4 = A_4 - Dot(A_4,U_1)*U_1 - Dot(A_4,U_2)*U_2 - Dot(A_4,U_3)*U_3\n",
    "U_4 = W_4/Magnitude(W_4)\n",
    "\n",
    "# Assemble the matrix Q\n",
    "\n",
    "Q = np.hstack((U_1,U_2,U_3,U_4))\n",
    "print(\"Q\")\n",
    "print(Q,'\\n')\n",
    "\n",
    "# Check that Q is orthogonal\n",
    "\n",
    "print(\"QTQ\")\n",
    "print(np.round(Q.transpose()@Q),'\\n')\n",
    "\n",
    "# Compute R\n",
    "\n",
    "R = Q.transpose()@A\n",
    "\n",
    "#  Check\n",
    "\n",
    "print(\"Q\")\n",
    "print(Q,'\\n')\n",
    "print(\"R\")\n",
    "print(np.round(R,8),'\\n')\n",
    "print(\"QR\")\n",
    "print(np.round(Q@R))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### QR factorization function\n",
    "\n",
    "It is helpful to write some code in a more general way so that we may easily carry out the factorization on an arbitrary matrix with linearly independent columns.  We will include in our function a simple check of the dimensions of the input matrix.  If the number of columns is greater than the number of rows, the columns cannot be linearly independent since the number of columns exceeds the dimension of the column space.  If this is the case, we will return an error message and not carry out any computations.  If the number of columns is not greater than the number of rows, we will leave it to the user of the function to check that the columns are linearly independent.  We can document this requirement with a comment in the code. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def QRFactorization(A):\n",
    "    # =============================================================================\n",
    "    # A is a Numpy array that represents a matrix of dimension m x n.\n",
    "    # QRFactorization returns matrices Q and R such that A=QR, Q is orthogonal\n",
    "    # and R is upper triangular.  The factorization is carried out using classical\n",
    "    # Gram-Schmidt and the results may suffer due to numerical instability.\n",
    "    # QRFactorization may not return correct results if the columns of A are \n",
    "    # linearly dependent.\n",
    "    # =============================================================================\n",
    "\n",
    "    # Check shape of A\n",
    "    if (A.shape[0] < A.shape[1]):\n",
    "        print(\"A must have more rows than columns for QR factorization.\")\n",
    "        return\n",
    "\n",
    "    m = A.shape[0]\n",
    "    n = A.shape[1]\n",
    "    \n",
    "    Q = np.zeros((m,n))\n",
    "    R = np.zeros((n,n))\n",
    "    \n",
    "    for i in range(n):\n",
    "        W = A[:,i:i+1]\n",
    "        for j in range(i):\n",
    "                W = W - Dot(A[:,i:i+1],Q[:,j:j+1])*Q[:,j:j+1]\n",
    "        Q[:,i:i+1] = W/Magnitude(W)\n",
    "        \n",
    "    R = Q.transpose()@A\n",
    "    \n",
    "    return (Q,R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One point to note here is that our function returns the pair of matrices $Q$ and $R$.  In order to receive both results, we need to assign a pair of names to the output.  This is exactly the same behavior as the SciPy function we used in Chapter 1 to compute the [PLU factorization](LU_Factorization.ipynb). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q,R = QRFactorization(A)\n",
    "print(Q)\n",
    "print('\\n')\n",
    "print(np.round(R,8))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 2:  Solving a linear system \n",
    "\n",
    "The orthogonalization behind the $QR$ factorization provides us another way to solve a linear system $AX=B$.  If we substitute $A=QR$, then multiply the equation by $Q^T$, we get $Q^TQRX = Q^TB$.  Once again $Q^TQ$ simplifies to $I$, so we are left with $RX=Q^TB$, which is a triangular system that can be solved easily by back substitution.\n",
    "\n",
    "Let's try it out on a $4\\times 4$ system.\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\left[ \\begin{array}{rrrr} \n",
    "2 & 3 & 0 & -1 \\\\ \n",
    "-1 & 0 & 2 & 0 \\\\ \n",
    "-1 & -1 & 4 & 2 \\\\\n",
    "0 & 3 & -3 & 2 \\\\\n",
    "\\end{array}\\right]X = \n",
    "\\left[ \\begin{array}{r} 0 \\\\ 1 \\\\ 2 \\\\ 5 \\end{array} \\right]\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from laguide import BackSubstitution\n",
    "\n",
    "A = np.array([[2, 3, 0, -1],[-1, 0, 2, 0],[-1, -1, 4, 2],[0, 3, -3, 2]])\n",
    "\n",
    "B = np.array([[0],[1],[2],[5]])\n",
    "X = np.array([[-1],[1],[0],[1]])\n",
    "Q,R = QRFactorization(A)\n",
    "C = Q.transpose()@B\n",
    "X = BackSubstitution(R,C)\n",
    "print(np.around(X,8),'\\n')\n",
    "\n",
    "## Check that AX = B\n",
    "print(np.around(A@X,8))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In some situations, we might find that we are solving several systems such as $AX=B_1$, $AX=B_2$, $AX=B_3$, ..., that involve the same matrix but different right hand sides.  In these situations it is useful to solve the systems with a factorization such as $QR$ or $LU$, because the factorization does not need to be recomputed for each system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### QR factorization with SciPy\n",
    "\n",
    "The SciPy function for finding $QR$ factorizations is called $\\texttt{qr}$.  The basic use is exactly the same as our $\\texttt{QRfactorization}$ function.  We apply $\\texttt{qr}$ to the matrix in the example and compare the results with our own computations.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.linalg as sla\n",
    "\n",
    "Qs,Rs = sla.qr(A)\n",
    "\n",
    "print(\"laguid Q\\n\",Q,'\\n',sep='')\n",
    "print(\"SciPy Q\\n\",Qs,'\\n',sep='')\n",
    "print(\"laguide R\\n\",np.around(R,8),'\\n',sep='')\n",
    "print(\"SciPy R\\n\",Rs,sep='')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although the results are displayed are close, we note that the entries of the SciPy are the negative of those we produced.  The SciPy function does not actually use the Gram Schmidt process, but instead approaches the factorization from an entirely different direction.    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "**Exercise 1:** \n",
    "\n",
    "($a$) Carry out the Gram-Schmidt algorithm on the following set of vectors to produce an orthonormal set $\\{U_1, U_2, U_3\\}$.  **Do not** use the $\\texttt{QRFactorization}$ function from this section or the SciPy $\\texttt{qr}$ function.  \n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "V_1 = \\left[ \\begin{array}{r} 0 \\\\ 2 \\\\ 0 \\\\ 1 \\end{array}\\right] \\hspace{0.7cm} \n",
    "V_2 = \\left[ \\begin{array}{r} 1 \\\\ 1 \\\\ 1 \\\\ 1 \\end{array}\\right] \\hspace{0.7cm}\n",
    "V_3 = \\left[ \\begin{array}{r} 2 \\\\ 2 \\\\ 0 \\\\ 1 \\end{array}\\right] \\hspace{0.7cm}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Code solution here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "($b$) Check your results by verifying that $\\{U_1, U_2, U_3\\}$ is an orthonormal set and that the span of $\\{U_1, U_2, U_3\\}$ equals the span of $\\{V_1, V_2, V_3\\}$.  (Reread [Linear_Combinations](Linear_Combinations.ipynb) if you are unsure how to verify the two spans are equal.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Code solution here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "($c$) Check your results against those obtained using the $\\texttt{QRFactorization}$ function from this section or the $\\texttt{qr}$ function from SciPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Code solution here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 2:** \n",
    "\n",
    "($a$) Predict what will happen if we attempt to find the QR factorization of matrix with *linearly dependent* columns.\n",
    "\n",
    "($b$) Try to compute the QR factorization on the following matrix with *linearly dependent* columns.  Try both the $\\texttt{QRFactorization}$ function from this section or the $\\texttt{qr}$ function from SciPy.\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "B = \\left[ \\begin{array}{rrr} \n",
    "1 & 3 & -1  \\\\ \n",
    "0 & -1 & 1  \\\\ \n",
    "2 & 2 & 2  \\\\\n",
    "1 & 1 & 1  \\\\\n",
    "1 & 0 & 2  \\\\\n",
    "\\end{array}\\right]\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Code solution here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 3:** If possible, find the $QR$ factorization of the matrix $A$. Try to find the matrices $Q$ and $R$ without using the $\\texttt{QRFactorization}$ function from this section, then check your result by verifying that $QR=A$. \n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "A = \\left[ \\begin{array}{rrr} \n",
    "1 & 3 & 0 & 2  \\\\ \n",
    "0 & 1 & 2 & 1 \\\\ \n",
    "2 & 1 & 2 & 1 \\\\\n",
    "1 & 0 & 1 & 3 \\\\\n",
    "\\end{array}\\right]\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Code solution here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 4:** Are all matrices which have $QR$ factorization invertible? Explain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Code solution here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 5:** Use the $QR$ factorization of $A$ to solve the given linear system $AX = B$, and verify the solution.\n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "AX = \\left[ \\begin{array}{rrrr} \n",
    "1 & 2 & 3  \\\\ \n",
    "0 & 3 & 2  \\\\ \n",
    "1 & 1 & 4  \\\\\n",
    "\\end{array}\\right]\n",
    "\\quad\\quad\n",
    "X = \n",
    "\\left[ \\begin{array}{r} 1 \\\\ 1 \\\\ 1 \\\\  \\end{array} \\right]\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Code solution here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 6:** $A$ and $B$ are two $n \\times n$ matrices. \n",
    "\n",
    "  ($a$) Given that the product $BA$ has $QR$ factorization, prove that the matrix $A$ also has a $QR$ factorization. \n",
    "\n",
    "\n",
    " ($b$) Can you think of matrices $A$ and $B$ such that $B$ has linearly dependent columns  but the product $BA$ has $QR$ factorization? Explain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Code solution here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 7:** If $A$ is an $n \\times n$ invertible matrix, explain why $A^TA$ is also invertible. "
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}