MGS_VS_HHT.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modified Gram-Schmidt (MGS) vs Householder’s Triangularization (HHT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison of Modified Gram-Schmidt (MGS) vs Householder’s Triangularization (HHT) in a computer assuming each computation is calculated on a computer which rounds all computed results to five digits of relative accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from typing import Tuple, Any\n",
    "from copy import deepcopy\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "A = np.array(\n",
    "    [\n",
    "        [0.70000, 0.70711],\n",
    "        [0.70001, 0.70711]\n",
    "    ]\n",
    ")\n",
    "\n",
    "I = np.array(\n",
    "    [\n",
    "        [1, 0],\n",
    "        [0, 1]\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def round_to_5(x : Any) -> Any:\n",
    "    return np.round(x, 5)\n",
    "\n",
    "def MGS(A : np.array) -> Tuple[np.array, np.array]:\n",
    "    \"\"\"\n",
    "    Returns Q, R for matrix A.\n",
    "    \"\"\"\n",
    "    m,n = A.shape\n",
    "    Q = np.zeros((m,n))\n",
    "    R = np.zeros((n,n))\n",
    "    V = np.zeros((m,n))\n",
    "\n",
    "    for i in range(n):\n",
    "        V[:,i] = deepcopy(A[:,i])\n",
    "    for i in range(n):\n",
    "        R[i,i] = round_to_5(np.linalg.norm(V[:,i],2)) \n",
    "        Q[:,i] = round_to_5(V[:,i]/R[i,i])\n",
    "\n",
    "        for j in range(i+1,n):\n",
    "            R[i,j] = round_to_5(np.dot(Q[:,i],V[:,j]))\n",
    "            V[:,j] = round_to_5(V[:,j] - R[i,j]*Q[:,i])\n",
    "\n",
    "    return Q, R\n",
    "\n",
    "def Householder(A: np.array) -> Tuple[np.array, np.array]:\n",
    "    \"\"\"\n",
    "    Takes in A and b to return Qtb and R using Householder Triangularization (Reduced).\n",
    "    \"\"\"\n",
    "    m,n = A.shape\n",
    "    R = deepcopy(A)\n",
    "    Q = np.eye(m)\n",
    "    for k in range(n):\n",
    "        x = deepcopy(R[k:, k])\n",
    "        e_1 = np.zeros_like(x); e_1[0] = 1\n",
    "        v_k = round_to_5(np.sign(x[0])* round_to_5(np.linalg.norm(x,2)) * e_1 + x)\n",
    "        v_k = np.reshape(v_k, (m-k,1))\n",
    "        norm_v_k = round_to_5(np.linalg.norm(v_k,2))\n",
    "        v_k = round_to_5(v_k / norm_v_k)\n",
    "        \n",
    "        R[k:, k:] = round_to_5(R[k:, k:] -  2 * v_k @ v_k.T @ R[k:, k:])\n",
    "        Q[:, k:] = round_to_5(Q[:, k:] - 2 * np.outer(np.dot(Q[:, k:], v_k), v_k))\n",
    "\n",
    "     # Trim Q to be reduced size (m x n)\n",
    "    Q_reduced = Q[:, :n]\n",
    "    \n",
    "    return Q_reduced, R[:n, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative condition number of A = 280016.2781264799\n",
      "=======================================================================\n",
      "Q matrix calculated by MGS :\n",
      " [[0.7071  1.     ]\n",
      " [0.70711 0.     ]]\n",
      "R matrix calculated by MGS :\n",
      " [[9.8996e-01 1.0000e+00]\n",
      " [0.0000e+00 1.0000e-05]]\n",
      "Loss of orthogonality ||I - Q_MGS.T Q_MGS|| = 0.9999904099667057\n",
      "=======================================================================\n",
      "\n",
      "=======================================================================\n",
      "Q matrix calculated by HHT :\n",
      " [[-0.70711  0.70712]\n",
      " [-0.70712 -0.7071 ]]\n",
      "Q matrix calculated by HHT :\n",
      " [[-9.89970e-01 -1.00001e+00]\n",
      " [-1.00000e-05  2.00000e-05]]\n",
      "Loss of orthogonality ||I - Q_HHT.T Q_HHT|| = 2.6894118325280188e-05\n",
      "=======================================================================\n",
      "\n"
     ]
    }
   ],
   "source": [
    "Q_MGS, R_MGS = MGS(A)\n",
    "Q_HHT, R_HHT = Householder(A)\n",
    "orthogonality_loss_MGS = np.linalg.norm(I - Q_MGS.T @ Q_MGS)\n",
    "orthogonality_loss_HHT = np.linalg.norm(I - Q_HHT.T @ Q_HHT)\n",
    "\n",
    "k_A = np.linalg.norm(A, 2)*np.linalg.norm(np.linalg.inv(A), 2)\n",
    "\n",
    "print(f\"Relative condition number of A = {k_A}\")\n",
    "print(\"=======================================================================\")\n",
    "print(f\"Q matrix calculated by MGS :\\n {Q_MGS}\")\n",
    "print(f\"R matrix calculated by MGS :\\n {R_MGS}\")\n",
    "print(f\"Loss of orthogonality ||I - Q_MGS.T Q_MGS|| = {orthogonality_loss_MGS}\")\n",
    "print(\"=======================================================================\\n\")\n",
    "\n",
    "print(\"=======================================================================\")\n",
    "print(f\"Q matrix calculated by HHT :\\n {Q_HHT}\")\n",
    "print(f\"Q matrix calculated by HHT :\\n {R_HHT}\")\n",
    "print(f\"Loss of orthogonality ||I - Q_HHT.T Q_HHT|| = {orthogonality_loss_HHT}\")\n",
    "print(\"=======================================================================\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see from the results above that Q_MGS computed by MGS has very high loss of orthogonality (0.99999) compared to Q_HHT ( loss =  2.68941e-05.). Therefore Householder algorithm is preferred on low precision machines.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}