{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# QR algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classical Gram-Schmidt algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qr(A):\n",
    "    n = A.shape[0]\n",
    "    Q = np.empty_like(A)\n",
    "    R = np.zeros((n,n))\n",
    "    for j in range(n):\n",
    "        Q[:,j] = A[:,j]\n",
    "        for i in range(j):\n",
    "            R[i,j] = Q[:,i].dot(A[:,j])\n",
    "            Q[:,j] -= R[i,j]*Q[:,i]\n",
    "        R[j,j] = np.linalg.norm(Q[:,j])\n",
    "        Q[:,j] /= R[j,j]\n",
    "    return Q,R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We generate a random matrix and test this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =\n",
      " [[0.65740643 0.80644261 0.58714836 0.5815085 ]\n",
      " [0.05973184 0.15401564 0.97914101 0.03653573]\n",
      " [0.6902891  0.9879373  0.73121611 0.73777795]\n",
      " [0.24723405 0.59180597 0.40238936 0.10926068]]\n",
      "Q =\n",
      " [[ 0.66633695 -0.41784852  0.05290702 -0.61530361]\n",
      " [ 0.06054327  0.24979359  0.9661765  -0.02099151]\n",
      " [ 0.69966631  0.06432462 -0.04504459  0.71014111]\n",
      " [ 0.2505926   0.87113037 -0.24834412 -0.3415559 ]]\n",
      "R =\n",
      " [[ 0.9865976   1.38621575  1.0629621   0.93327088]\n",
      " [ 0.          0.28059011  0.39681286 -0.09121849]\n",
      " [ 0.          0.          0.84421894  0.00569869]\n",
      " [ 0.          0.          0.          0.1280366 ]]\n",
      "A-QR =\n",
      " [[0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
      " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
      " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
      " [0.00000000e+00 0.00000000e+00 5.55111512e-17 0.00000000e+00]]\n"
     ]
    }
   ],
   "source": [
    "n = 4\n",
    "A = np.random.rand(n,n)\n",
    "print('A =\\n',A)\n",
    "Q,R = qr(A)\n",
    "print('Q =\\n',Q)\n",
    "print('R =\\n',R)\n",
    "print('A-QR =\\n',A-Q.dot(R))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modified Gram-Schmidt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mqr(A):\n",
    "    n = A.shape[0]\n",
    "    Q = A.copy()\n",
    "    R = np.zeros((n,n))\n",
    "    for i in range(n):\n",
    "        R[i,i] = np.linalg.norm(Q[:,i])\n",
    "        Q[:,i] /= R[i,i]\n",
    "        for j in range(i+1,n):\n",
    "            R[i,j] = Q[:,i].dot(Q[:,j])\n",
    "            Q[:,j] -= R[i,j]*Q[:,i]\n",
    "    return Q,R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Q =\n",
      " [[ 0.66633695 -0.41784852  0.05290702 -0.61530361]\n",
      " [ 0.06054327  0.24979359  0.9661765  -0.02099151]\n",
      " [ 0.69966631  0.06432462 -0.04504459  0.71014111]\n",
      " [ 0.2505926   0.87113037 -0.24834412 -0.3415559 ]]\n",
      "R =\n",
      " [[ 0.9865976   1.38621575  1.0629621   0.93327088]\n",
      " [ 0.          0.28059011  0.39681286 -0.09121849]\n",
      " [ 0.          0.          0.84421894  0.00569869]\n",
      " [ 0.          0.          0.          0.1280366 ]]\n",
      "A-QR =\n",
      " [[0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
      " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
      " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
      " [0.00000000e+00 0.00000000e+00 5.55111512e-17 0.00000000e+00]]\n"
     ]
    }
   ],
   "source": [
    "Q,R = mqr(A)\n",
    "print('Q =\\n',Q)\n",
    "print('R =\\n',R)\n",
    "print('A-QR =\\n',A-Q.dot(R))"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}