tspspi/qrdecomposition.ipynb

## qrdecomposition.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3aae8f85",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bc603ac",
   "metadata": {},
   "source": [
    "# QR decomposition using Givens rotations\n",
    "\n",
    "In this notebook there will be really simple demonstration on how to perform QR decomposition (as described in https://www.tspi.at/2021/12/08/qrgivens.html) using Givens rotations. Note that this also already shows how it's affected by rounding precision."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b78c67dc",
   "metadata": {},
   "source": [
    "## A simple matrix class\n",
    "\n",
    "First lets define a really simple and basic matrix class - this is the worst possible implementation when it comes to performance though but it suites well for demonstration purposes. It just implements the most basic methods:\n",
    "\n",
    "* Creation (and list based initialization)\n",
    "* Getting and setting values without knowing anything about internal storage\n",
    "* Copying\n",
    "* Getting name and description\n",
    "* Transposing\n",
    "* Adding\n",
    "* Subtracting\n",
    "* Multiplying with other matrices (including vectors) and scalars\n",
    "\n",
    "All indexing is done in the row-column order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9a9935d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Matrix:\n",
    "    def __init__(self, rows, cols, *, rounding = 8, initList = None, name = None, description = None):\n",
    "        self.name = name\n",
    "        self.description = description\n",
    "        self.rows = rows\n",
    "        self.cols = cols\n",
    "        self.data = [ 0 ] * rows\n",
    "        self.rounding = rounding\n",
    "        for row in range(rows):\n",
    "            self.data[row] = [ 0.0 ] * cols\n",
    "        # The following would NOT work since we create the same copy of the\n",
    "        # lists and reference them\n",
    "        # self.data = [ [ 0.0 ] * cols ] * rows\n",
    "        if initList:\n",
    "            if not isinstance(initList, list):\n",
    "                raise ValueError(\"List initializer has to be a list ...\")\n",
    "            # Two ways: Either one list with rows * cols entries just ordered in row major\n",
    "            # order (first row, second row, etc.) or nested list\n",
    "            if len(initList) == (rows * cols):\n",
    "                for row in range(rows):\n",
    "                    for col in range(cols):\n",
    "                        try:\n",
    "                            k = float(initList[col + row * cols])\n",
    "                            self.data[row][col] = round(k, self.rounding)\n",
    "                        except ValueError:\n",
    "                            raise ValueError(\"Initializer list has to contain numbers\")\n",
    "            else:\n",
    "                if len(initList) != rows:\n",
    "                    raise ValueError(\"Row mismatch in initializer list\")\n",
    "                for row in range(rows):\n",
    "                    if not isinstance(initList[row], list):\n",
    "                        raise ValueError(\"Columns have to be lists in initializer\")\n",
    "                    if len(initList[row]) != cols:\n",
    "                        raise ValueError(\"Column mismatch in initializer list\")\n",
    "                    for col in range(cols):\n",
    "                        try:\n",
    "                            k = float(initList[row][col])\n",
    "                            self.data[row][col] = round(k, self.rounding)\n",
    "                        except ValueError:\n",
    "                            raise ValueError(\"Initializer list has to contain numbers\")\n",
    "    def get(self, row, col):\n",
    "        if (row < 0) or (row >= self.rows):\n",
    "            raise IndexError(\"Row {} invalid (matrix is {} x {})\".format(row, self.rows, self.cols))\n",
    "        if (col < 0) or (col >= self.cols):\n",
    "            raise IndexError(\"Column {} invalid (matrix is {} x {})\".format(col, self.rows, self.cols))\n",
    "        return self.data[row][col]\n",
    "    def set(self, row, col, value):\n",
    "        if (row < 0) or (row >= self.rows):\n",
    "            raise IndexError(\"Row {} invalid (matrix is {} x {})\".format(row, self.rows, self.cols))\n",
    "        if (col < 0) or (col >= self.cols):\n",
    "            raise IndexError(\"Column {} invalid (matrix is {} x {})\".format(col, self.rows, self.cols))\n",
    "        self.data[row][col] = round(value, self.rounding)\n",
    "    def copy(self, *, name = None, description = None):\n",
    "        result = Matrix(self.rows, self.cols, name = name, description = description)\n",
    "        for row in range(self.rows):\n",
    "            for col in range(self.cols):\n",
    "                result.set(row, col, self.get(row, col))\n",
    "        return result\n",
    "    def getName(self):\n",
    "        return self.name\n",
    "    def getDescription(self):\n",
    "        return self.description\n",
    "    def __str__(self):\n",
    "        resString = \"\"\n",
    "        if self.name:\n",
    "            resString = resString + \"{} = \".format(self.name)\n",
    "        for row in range(self.rows):\n",
    "            for col in range(self.cols):\n",
    "                resString = resString + \"{}\".format(self.data[row][col]) + \"\\t\"\n",
    "            resString = resString + \"\\n\"\n",
    "        return resString\n",
    "    def transpose(self, *, name = None, description = None):\n",
    "        result = Matrix(self.cols, self.rows, name = name, description = description)\n",
    "        for col in range(self.cols):\n",
    "            for row in range(self.rows):\n",
    "                result.set(col, row, self.get(row, col))\n",
    "        return result\n",
    "    def add(self, other, *, name = None, description = None):\n",
    "        if not isinstance(other, Matrix):\n",
    "            raise ValueError(\"Adding matrix with non matrix is not supported\")\n",
    "        if (other.cols != self.cols) or (other.rows != self.rows):\n",
    "            raise ValueError(\"Addition requires same dimensions\")\n",
    "        result = Matrix(self.rows, self.cols, name = name, description = description)\n",
    "        for col in range(self.cols):\n",
    "            for row in range(self.rows):\n",
    "                result.set(row, col, round(self.get(row, col) + other.get(row, col), self.rounding))\n",
    "        return result\n",
    "    def sub(self, other, *, name = None, description = None):\n",
    "        if not isinstance(other, Matrix):\n",
    "            raise ValueError(\"Adding matrix with non matrix is not supported\")\n",
    "        if (other.cols != self.cols) or (other.rows != self.rows):\n",
    "            raise ValueError(\"Addition requires same dimensions\")\n",
    "        result = Matrix(self.rows, self.cols, name = name, description = description)\n",
    "        for col in range(self.cols):\n",
    "            for row in range(self.rows):\n",
    "                result.set(row, col, round(self.get(row, col) - other.get(row, col), self.rounding))\n",
    "        return result\n",
    "    def multiply(self, other, *, name = None, description = None):\n",
    "        # Calculate self * other\n",
    "        if isinstance(other, Matrix):\n",
    "            if not self.cols == other.rows:\n",
    "                raise ValueError(\"Colums of left side {} have to match rows of right side {}\".format(self.cols, other.rows))\n",
    "            result = Matrix(self.rows, other.cols, name = name, description = description)\n",
    "            for row in range(self.rows):\n",
    "                for col in range(other.cols):\n",
    "                    t = 0\n",
    "                    for k in range(self.cols):\n",
    "                        t = t + self.get(row, k) * other.get(k, col)\n",
    "                    t = round(t, 4)\n",
    "                    result.set(row, col, round(t, self.rounding))\n",
    "            return result\n",
    "        else:\n",
    "            try:\n",
    "                otherFloat = float(other)\n",
    "                result = Matrix(self.rows, self.cols, name = name, description = description)\n",
    "                for r in range(self.rows):\n",
    "                    for c in range(self.cols):\n",
    "                        result.set(r,c, round(self.get(r,c,)*otherFloat, self.rounding))\n",
    "                return result\n",
    "            except ValueError:\n",
    "                raise ValueError(\"Matrix has to be multiplied by matrix or scalar\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f32d5200",
   "metadata": {},
   "source": [
    "As a helper let's also define a simple identity matrix class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "654a576a",
   "metadata": {},
   "outputs": [],
   "source": [
    "class IdentityMatrix(Matrix):\n",
    "    def __init__(self, size, *, name = None, description = None):\n",
    "        super().__init__(size, size, name = name, description = description)\n",
    "        for k in range(size):\n",
    "            self.set(k,k,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f692886a",
   "metadata": {},
   "source": [
    "## The QR decomposition\n",
    "\n",
    "To perform QR decomposition using Givens rotations let's first define a class that creates the $G_{i,j}$ matrix when we also tell it the values of $x_i$ and $x_k$ in our given column where $i$ and $k$ are the row indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9834c499",
   "metadata": {},
   "outputs": [],
   "source": [
    "class GivensRotationMatrix(IdentityMatrix):\n",
    "    def __init__(self, size, i, j, x1, x2, *, name = None, description = None):\n",
    "        super().__init__(size, name = name, description = description)\n",
    "        if (i >= size) or (i < 0) or (j >= size) or (j < 0):\n",
    "            raise ValueError(\"Invalid first index i or j\")\n",
    "        # Calculate s and c\n",
    "        s = 1.0 / math.sqrt(x1 * x1 / (x2 * x2) + 1)\n",
    "        c = s * x1 / x2\n",
    "        self.set(i,i,round(c, self.rounding))\n",
    "        self.set(j,i,round(-1 * s, self.rounding))\n",
    "        self.set(i,j,round(s, self.rounding))\n",
    "        self.set(j,j,round(c, self.rounding))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b96346c6",
   "metadata": {},
   "source": [
    "Now we know that we should apply the transformations to zero out the lower triangular area starting from the lowest row in the leftmost column proceeding upwards into the first subdiagonal, then continuing the same process from bottom to top in the second column, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6b90b208",
   "metadata": {},
   "outputs": [],
   "source": [
    "def QRDecomposeWithGivens(matIn, *, traceOutput = False):\n",
    "    matR = matIn.copy()\n",
    "    matQ = IdentityMatrix(matIn.rows)\n",
    "    if traceOutput:\n",
    "        print(\"Input into QR decomposition:\")\n",
    "        print(matIn)\n",
    "    for col in range(matIn.cols):\n",
    "        for row in range(matIn.rows-1, col, -1):\n",
    "            if traceOutput:\n",
    "                x1 = matR.get(row-1, col)\n",
    "                x2 = matR.get(row,col)\n",
    "            rotationMatrix = GivensRotationMatrix(matIn.rows, row-1, row, matR.get(row-1, col), matR.get(row,col))\n",
    "            matR = rotationMatrix.multiply(matR)\n",
    "            matQ = matQ.multiply(rotationMatrix.transpose())\n",
    "            if traceOutput:\n",
    "                print(\"--- Next step ---\")\n",
    "                print(\"Rotation matrix for col {}, row {}; x{}: {}, x{}: {}\".format(col, row, row-1, x1, row, x2))\n",
    "                print(rotationMatrix)\n",
    "                print(\"New R\")\n",
    "                print(matR)\n",
    "                print(\"New Q\")\n",
    "                print(matQ)\n",
    "\n",
    "    return (matQ, matR)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08b1cfff",
   "metadata": {},
   "source": [
    "## A simple test\n",
    "\n",
    "To verify the decomposition works let's use a simple test using a matrix that has also been found in some textbook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "cda3ed6a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8147\t0.0975\t0.1576\t\n",
      "0.9058\t0.2785\t0.9706\t\n",
      "0.127\t0.5469\t0.9572\t\n",
      "0.9134\t0.9575\t0.4854\t\n",
      "0.6324\t0.9649\t0.8003\t\n",
      "\n"
     ]
    }
   ],
   "source": [
    "testMatrix1 = Matrix(5, 3, initList = [ [ 0.8147, 0.0975, 0.1576 ], [ 0.9058, 0.2785, 0.9706 ], [ 0.1270, 0.5469, 0.9572 ], [ 0.9134, 0.9575, 0.4854 ], [ 0.6324, 0.9649, 0.8003 ] ])\n",
    "print(testMatrix1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01628469",
   "metadata": {},
   "source": [
    "Now apply our QR decomposition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "72a8e3da",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Q:\n",
      "0.4927\t-0.4807\t-0.1779\t-0.7032\t0.0\t\n",
      "0.5477\t-0.3583\t0.5778\t0.4824\t-0.0706\t\n",
      "0.0768\t0.4754\t0.6342\t-0.4316\t0.4235\t\n",
      "0.5524\t0.3391\t-0.481\t0.2769\t0.5216\t\n",
      "0.3824\t0.5473\t-0.0312\t-0.0984\t-0.7373\t\n",
      "\n",
      "R:\n",
      "1.6536\t1.1405\t1.257\t\n",
      "0.0\t0.9662\t0.6341\t\n",
      "0.0\t0.0\t0.8815\t\n",
      "0.0\t0.0\t-0.0\t\n",
      "0.0\t0.0\t0.0\t\n",
      "\n"
     ]
    }
   ],
   "source": [
    "q,r = QRDecomposeWithGivens(testMatrix1, traceOutput = False)\n",
    "\n",
    "print(\"Q:\")\n",
    "print(q)\n",
    "print(\"R:\")\n",
    "print(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "107c3d02",
   "metadata": {},
   "source": [
    "And calculate the product $Q*R$ again that should yield our input (up to rounding errors). In the end also calculate the difference to make comparison easier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6e07cb68",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8147\t0.0975\t0.1577\t\n",
      "0.9057\t0.2785\t0.9706\t\n",
      "0.127\t0.5469\t0.957\t\n",
      "0.9134\t0.9577\t0.4854\t\n",
      "0.6323\t0.9649\t0.8002\t\n",
      "\n"
     ]
    }
   ],
   "source": [
    "qr = q.multiply(r)\n",
    "print(qr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "760a34dd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\t0.0\t-0.0001\t\n",
      "0.0001\t0.0\t0.0\t\n",
      "0.0\t0.0\t0.0002\t\n",
      "0.0\t-0.0002\t0.0\t\n",
      "0.0001\t0.0\t0.0001\t\n",
      "\n"
     ]
    }
   ],
   "source": [
    "difference = testMatrix1.sub(qr)\n",
    "print(difference)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9c91011",
   "metadata": {},
   "source": [
    "As one can see multiplying the Q and R matrices again yields the initial matrix. To verify if we have an (up to rounding errors) orthogonal matrix let's check if $Q^T * Q = \\mathfrak{1}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "5dcb4570",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\t0.0\t-0.0001\t-0.0001\t0.0\t\n",
      "0.0\t1.0\t-0.0002\t0.0\t-0.0\t\n",
      "-0.0001\t-0.0002\t1.0\t-0.0\t-0.0001\t\n",
      "-0.0001\t0.0\t-0.0\t0.9998\t0.0001\t\n",
      "0.0\t-0.0\t-0.0001\t0.0001\t1.0\t\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(q.transpose().multiply(q))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "555c1134",
   "metadata": {},
   "source": [
    "Up to rounding errors this yields acceptable results as an orthogonal matrix. Thus our decomposition seems to really work as expected."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec325a34",
   "metadata": {},
   "source": [
    "## Including trace output\n",
    "\n",
    "To make verification of implementations easier the same QR decomposition can trace all intermediate steps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "4591c08c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Input into QR decomposition:\n",
      "0.8147\t0.0975\t0.1576\t\n",
      "0.9058\t0.2785\t0.9706\t\n",
      "0.127\t0.5469\t0.9572\t\n",
      "0.9134\t0.9575\t0.4854\t\n",
      "0.6324\t0.9649\t0.8003\t\n",
      "\n",
      "--- Next step ---\n",
      "Rotation matrix for col 0, row 4; x3: 0.9134, x4: 0.6324\n",
      "1\t0.0\t0.0\t0.0\t0.0\t\n",
      "0.0\t1\t0.0\t0.0\t0.0\t\n",
      "0.0\t0.0\t1\t0.0\t0.0\t\n",
      "0.0\t0.0\t0.0\t0.82217275\t0.56923806\t\n",
      "0.0\t0.0\t0.0\t-0.56923806\t0.82217275\t\n",
      "\n",
      "New R\n",
      "0.8147\t0.0975\t0.1576\t\n",
      "0.9058\t0.2785\t0.9706\t\n",
      "0.127\t0.5469\t0.9572\t\n",
      "1.111\t1.3365\t0.8546\t\n",
      "0.0\t0.2483\t0.3817\t\n",
      "\n",
      "New Q\n",
      "1.0\t0.0\t0.0\t0.0\t0.0\t\n",
      "0.0\t1.0\t0.0\t0.0\t0.0\t\n",
      "0.0\t0.0\t1.0\t0.0\t0.0\t\n",
      "0.0\t0.0\t0.0\t0.8222\t-0.5692\t\n",
      "0.0\t0.0\t0.0\t0.5692\t0.8222\t\n",
      "\n",
      "--- Next step ---\n",
      "Rotation matrix for col 0, row 3; x2: 0.127, x3: 1.111\n",
      "1\t0.0\t0.0\t0.0\t0.0\t\n",
      "0.0\t1\t0.0\t0.0\t0.0\t\n",
      "0.0\t0.0\t0.11357181\t0.99352979\t0.0\t\n",
      "0.0\t0.0\t-0.99352979\t0.11357181\t0.0\t\n",
      "0.0\t0.0\t0.0\t0.0\t1\t\n",
      "\n",
      "New R\n",
      "0.8147\t0.0975\t0.1576\t\n",
      "0.9058\t0.2785\t0.9706\t\n",
      "1.1182\t1.39\t0.9578\t\n",
      "-0.0\t-0.3916\t-0.8539\t\n",
      "0.0\t0.2483\t0.3817\t\n",
      "\n",
      "New Q\n",
      "1.0\t0.0\t0.0\t0.0\t0.0\t\n",
      "0.0\t1.0\t0.0\t0.0\t0.0\t\n",
      "0.0\t0.0\t0.1136\t-0.9935\t0.0\t\n",
      "0.0\t0.0\t0.8169\t0.0934\t-0.5692\t\n",
      "0.0\t0.0\t0.5655\t0.0646\t0.8222\t\n",
      "\n",
      "--- Next step ---\n",
      "Rotation matrix for col 0, row 2; x1: 0.9058, x2: 1.1182\n",
      "1\t0.0\t0.0\t0.0\t0.0\t\n",
      "0.0\t0.62944608\t0.77704417\t0.0\t0.0\t\n",
      "0.0\t-0.77704417\t0.62944608\t0.0\t0.0\t\n",
      "0.0\t0.0\t0.0\t1\t0.0\t\n",
      "0.0\t0.0\t0.0\t0.0\t1\t\n",
      "\n",
      "New R\n",
      "0.8147\t0.0975\t0.1576\t\n",
      "1.439\t1.2554\t1.3552\t\n",
      "-0.0\t0.6585\t-0.1513\t\n",
      "0.0\t-0.3916\t-0.8539\t\n",
      "0.0\t0.2483\t0.3817\t\n",
      "\n",
      "New Q\n",
      "1.0\t0.0\t0.0\t0.0\t0.0\t\n",
      "0.0\t0.6294\t-0.777\t0.0\t0.0\t\n",
      "0.0\t0.0883\t0.0715\t-0.9935\t0.0\t\n",
      "0.0\t0.6348\t0.5142\t0.0934\t-0.5692\t\n",
      "0.0\t0.4394\t0.356\t0.0646\t0.8222\t\n",
      "\n",
      "--- Next step ---\n",
      "Rotation matrix for col 0, row 1; x0: 0.8147, x1: 1.439\n",
      "0.49267685\t0.87021234\t0.0\t0.0\t0.0\t\n",
      "-0.87021234\t0.49267685\t0.0\t0.0\t0.0\t\n",
      "0.0\t0.0\t1\t0.0\t0.0\t\n",
      "0.0\t0.0\t0.0\t1\t0.0\t\n",
      "0.0\t0.0\t0.0\t0.0\t1\t\n",
      "\n",
      "New R\n",
      "1.6536\t1.1405\t1.257\t\n",
      "-0.0\t0.5337\t0.5305\t\n",
      "0.0\t0.6585\t-0.1513\t\n",
      "0.0\t-0.3916\t-0.8539\t\n",
      "0.0\t0.2483\t0.3817\t\n",
      "\n",
      "New Q\n",
      "0.4927\t-0.8702\t0.0\t0.0\t0.0\t\n",
      "0.5477\t0.3101\t-0.777\t0.0\t0.0\t\n",
      "0.0768\t0.0435\t0.0715\t-0.9935\t0.0\t\n",
      "0.5524\t0.3128\t0.5142\t0.0934\t-0.5692\t\n",
      "0.3824\t0.2165\t0.356\t0.0646\t0.8222\t\n",
      "\n",
      "--- Next step ---\n",
      "Rotation matrix for col 1, row 4; x3: -0.3916, x4: 0.2483\n",
      "1\t0.0\t0.0\t0.0\t0.0\t\n",
      "0.0\t1\t0.0\t0.0\t0.0\t\n",
      "0.0\t0.0\t1\t0.0\t0.0\t\n",
      "0.0\t0.0\t0.0\t-0.8445395\t0.53549326\t\n",
      "0.0\t0.0\t0.0\t-0.53549326\t-0.8445395\t\n",
      "\n",
      "New R\n",
      "1.6536\t1.1405\t1.257\t\n",
      "0.0\t0.5337\t0.5305\t\n",
      "0.0\t0.6585\t-0.1513\t\n",
      "0.0\t0.4637\t0.9256\t\n",
      "0.0\t0.0\t0.1349\t\n",
      "\n",
      "New Q\n",
      "0.4927\t-0.8702\t0.0\t0.0\t0.0\t\n",
      "0.5477\t0.3101\t-0.777\t0.0\t0.0\t\n",
      "0.0768\t0.0435\t0.0715\t0.839\t0.532\t\n",
      "0.5524\t0.3128\t0.5142\t-0.3837\t0.4307\t\n",
      "0.3824\t0.2165\t0.356\t0.3857\t-0.729\t\n",
      "\n",
      "--- Next step ---\n",
      "Rotation matrix for col 1, row 3; x2: 0.6585, x3: 0.4637\n",
      "1\t0.0\t0.0\t0.0\t0.0\t\n",
      "0.0\t1\t0.0\t0.0\t0.0\t\n",
      "0.0\t0.0\t0.81762457\t0.57575173\t0.0\t\n",
      "0.0\t0.0\t-0.57575173\t0.81762457\t0.0\t\n",
      "0.0\t0.0\t0.0\t0.0\t1\t\n",
      "\n",
      "New R\n",
      "1.6536\t1.1405\t1.257\t\n",
      "0.0\t0.5337\t0.5305\t\n",
      "0.0\t0.8054\t0.4092\t\n",
      "0.0\t-0.0\t0.8439\t\n",
      "0.0\t0.0\t0.1349\t\n",
      "\n",
      "New Q\n",
      "0.4927\t-0.8702\t0.0\t0.0\t0.0\t\n",
      "0.5477\t0.3101\t-0.6353\t0.4474\t0.0\t\n",
      "0.0768\t0.0435\t0.5415\t0.6448\t0.532\t\n",
      "0.5524\t0.3128\t0.1995\t-0.6098\t0.4307\t\n",
      "0.3824\t0.2165\t0.5131\t0.1104\t-0.729\t\n",
      "\n",
      "--- Next step ---\n",
      "Rotation matrix for col 1, row 2; x1: 0.5337, x2: 0.8054\n",
      "1\t0.0\t0.0\t0.0\t0.0\t\n",
      "0.0\t0.55238123\t0.83359161\t0.0\t0.0\t\n",
      "0.0\t-0.83359161\t0.55238123\t0.0\t0.0\t\n",
      "0.0\t0.0\t0.0\t1\t0.0\t\n",
      "0.0\t0.0\t0.0\t0.0\t1\t\n",
      "\n",
      "New R\n",
      "1.6536\t1.1405\t1.257\t\n",
      "0.0\t0.9662\t0.6341\t\n",
      "0.0\t0.0\t-0.2162\t\n",
      "0.0\t0.0\t0.8439\t\n",
      "0.0\t0.0\t0.1349\t\n",
      "\n",
      "New Q\n",
      "0.4927\t-0.4807\t0.7254\t0.0\t0.0\t\n",
      "0.5477\t-0.3583\t-0.6094\t0.4474\t0.0\t\n",
      "0.0768\t0.4754\t0.2629\t0.6448\t0.532\t\n",
      "0.5524\t0.3391\t-0.1505\t-0.6098\t0.4307\t\n",
      "0.3824\t0.5473\t0.103\t0.1104\t-0.729\t\n",
      "\n",
      "--- Next step ---\n",
      "Rotation matrix for col 2, row 4; x3: 0.8439, x4: 0.1349\n",
      "1\t0.0\t0.0\t0.0\t0.0\t\n",
      "0.0\t1\t0.0\t0.0\t0.0\t\n",
      "0.0\t0.0\t1\t0.0\t0.0\t\n",
      "0.0\t0.0\t0.0\t0.98746326\t0.15784903\t\n",
      "0.0\t0.0\t0.0\t-0.15784903\t0.98746326\t\n",
      "\n",
      "New R\n",
      "1.6536\t1.1405\t1.257\t\n",
      "0.0\t0.9662\t0.6341\t\n",
      "0.0\t0.0\t-0.2162\t\n",
      "0.0\t0.0\t0.8546\t\n",
      "0.0\t0.0\t-0.0\t\n",
      "\n",
      "New Q\n",
      "0.4927\t-0.4807\t0.7254\t0.0\t0.0\t\n",
      "0.5477\t-0.3583\t-0.6094\t0.4418\t-0.0706\t\n",
      "0.0768\t0.4754\t0.2629\t0.7207\t0.4235\t\n",
      "0.5524\t0.3391\t-0.1505\t-0.5342\t0.5216\t\n",
      "0.3824\t0.5473\t0.103\t-0.0061\t-0.7373\t\n",
      "\n",
      "--- Next step ---\n",
      "Rotation matrix for col 2, row 3; x2: -0.2162, x3: 0.8546\n",
      "1\t0.0\t0.0\t0.0\t0.0\t\n",
      "0.0\t1\t0.0\t0.0\t0.0\t\n",
      "0.0\t0.0\t-0.24525723\t0.96945804\t0.0\t\n",
      "0.0\t0.0\t-0.96945804\t-0.24525723\t0.0\t\n",
      "0.0\t0.0\t0.0\t0.0\t1\t\n",
      "\n",
      "New R\n",
      "1.6536\t1.1405\t1.257\t\n",
      "0.0\t0.9662\t0.6341\t\n",
      "0.0\t0.0\t0.8815\t\n",
      "0.0\t0.0\t-0.0\t\n",
      "0.0\t0.0\t0.0\t\n",
      "\n",
      "New Q\n",
      "0.4927\t-0.4807\t-0.1779\t-0.7032\t0.0\t\n",
      "0.5477\t-0.3583\t0.5778\t0.4824\t-0.0706\t\n",
      "0.0768\t0.4754\t0.6342\t-0.4316\t0.4235\t\n",
      "0.5524\t0.3391\t-0.481\t0.2769\t0.5216\t\n",
      "0.3824\t0.5473\t-0.0312\t-0.0984\t-0.7373\t\n",
      "\n"
     ]
    }
   ],
   "source": [
    "q,r = QRDecomposeWithGivens(testMatrix1, traceOutput = True)"
   ]
  }
 ],
 "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
	{
	"cells": [
	{
	"cell_type": "code",
	"execution_count": 1,
	"id": "3aae8f85",
	"metadata": {},
	"outputs": [],
	"source": [
	"import math"
	]
	},
	{
	"cell_type": "markdown",
	"id": "9bc603ac",
	"metadata": {},
	"source": [
	"# QR decomposition using Givens rotations\n",
	"\n",
	"In this notebook there will be really simple demonstration on how to perform QR decomposition (as described in https://www.tspi.at/2021/12/08/qrgivens.html) using Givens rotations. Note that this also already shows how it's affected by rounding precision."
	]
	},
	{
	"cell_type": "markdown",
	"id": "b78c67dc",
	"metadata": {},
	"source": [
	"## A simple matrix class\n",
	"\n",
	"First lets define a really simple and basic matrix class - this is the worst possible implementation when it comes to performance though but it suites well for demonstration purposes. It just implements the most basic methods:\n",
	"\n",
	"* Creation (and list based initialization)\n",
	"* Getting and setting values without knowing anything about internal storage\n",
	"* Copying\n",
	"* Getting name and description\n",
	"* Transposing\n",
	"* Adding\n",
	"* Subtracting\n",
	"* Multiplying with other matrices (including vectors) and scalars\n",
	"\n",
	"All indexing is done in the row-column order."
	]
	},
	{
	"cell_type": "code",
	"execution_count": 2,
	"id": "9a9935d0",
	"metadata": {},
	"outputs": [],
	"source": [
	"class Matrix:\n",
	" def __init__(self, rows, cols, *, rounding = 8, initList = None, name = None, description = None):\n",
	" self.name = name\n",
	" self.description = description\n",
	" self.rows = rows\n",
	" self.cols = cols\n",
	" self.data = [ 0 ] * rows\n",
	" self.rounding = rounding\n",
	" for row in range(rows):\n",
	" self.data[row] = [ 0.0 ] * cols\n",
	" # The following would NOT work since we create the same copy of the\n",
	" # lists and reference them\n",
	" # self.data = [ [ 0.0 ] * cols ] * rows\n",
	" if initList:\n",
	" if not isinstance(initList, list):\n",
	" raise ValueError(\"List initializer has to be a list ...\")\n",
	" # Two ways: Either one list with rows * cols entries just ordered in row major\n",
	" # order (first row, second row, etc.) or nested list\n",
	" if len(initList) == (rows * cols):\n",
	" for row in range(rows):\n",
	" for col in range(cols):\n",
	" try:\n",
	" k = float(initList[col + row * cols])\n",
	" self.data[row][col] = round(k, self.rounding)\n",
	" except ValueError:\n",
	" raise ValueError(\"Initializer list has to contain numbers\")\n",
	" else:\n",
	" if len(initList) != rows:\n",
	" raise ValueError(\"Row mismatch in initializer list\")\n",
	" for row in range(rows):\n",
	" if not isinstance(initList[row], list):\n",
	" raise ValueError(\"Columns have to be lists in initializer\")\n",
	" if len(initList[row]) != cols:\n",
	" raise ValueError(\"Column mismatch in initializer list\")\n",
	" for col in range(cols):\n",
	" try:\n",
	" k = float(initList[row][col])\n",
	" self.data[row][col] = round(k, self.rounding)\n",
	" except ValueError:\n",
	" raise ValueError(\"Initializer list has to contain numbers\")\n",
	" def get(self, row, col):\n",
	" if (row < 0) or (row >= self.rows):\n",
	" raise IndexError(\"Row {} invalid (matrix is {} x {})\".format(row, self.rows, self.cols))\n",
	" if (col < 0) or (col >= self.cols):\n",
	" raise IndexError(\"Column {} invalid (matrix is {} x {})\".format(col, self.rows, self.cols))\n",
	" return self.data[row][col]\n",
	" def set(self, row, col, value):\n",
	" if (row < 0) or (row >= self.rows):\n",
	" raise IndexError(\"Row {} invalid (matrix is {} x {})\".format(row, self.rows, self.cols))\n",
	" if (col < 0) or (col >= self.cols):\n",
	" raise IndexError(\"Column {} invalid (matrix is {} x {})\".format(col, self.rows, self.cols))\n",
	" self.data[row][col] = round(value, self.rounding)\n",
	" def copy(self, *, name = None, description = None):\n",
	" result = Matrix(self.rows, self.cols, name = name, description = description)\n",
	" for row in range(self.rows):\n",
	" for col in range(self.cols):\n",
	" result.set(row, col, self.get(row, col))\n",
	" return result\n",
	" def getName(self):\n",
	" return self.name\n",
	" def getDescription(self):\n",
	" return self.description\n",
	" def __str__(self):\n",
	" resString = \"\"\n",
	" if self.name:\n",
	" resString = resString + \"{} = \".format(self.name)\n",
	" for row in range(self.rows):\n",
	" for col in range(self.cols):\n",
	" resString = resString + \"{}\".format(self.data[row][col]) + \"\\t\"\n",
	" resString = resString + \"\\n\"\n",
	" return resString\n",
	" def transpose(self, *, name = None, description = None):\n",
	" result = Matrix(self.cols, self.rows, name = name, description = description)\n",
	" for col in range(self.cols):\n",
	" for row in range(self.rows):\n",
	" result.set(col, row, self.get(row, col))\n",
	" return result\n",
	" def add(self, other, *, name = None, description = None):\n",
	" if not isinstance(other, Matrix):\n",
	" raise ValueError(\"Adding matrix with non matrix is not supported\")\n",
	" if (other.cols != self.cols) or (other.rows != self.rows):\n",
	" raise ValueError(\"Addition requires same dimensions\")\n",
	" result = Matrix(self.rows, self.cols, name = name, description = description)\n",
	" for col in range(self.cols):\n",
	" for row in range(self.rows):\n",
	" result.set(row, col, round(self.get(row, col) + other.get(row, col), self.rounding))\n",
	" return result\n",
	" def sub(self, other, *, name = None, description = None):\n",
	" if not isinstance(other, Matrix):\n",
	" raise ValueError(\"Adding matrix with non matrix is not supported\")\n",
	" if (other.cols != self.cols) or (other.rows != self.rows):\n",
	" raise ValueError(\"Addition requires same dimensions\")\n",
	" result = Matrix(self.rows, self.cols, name = name, description = description)\n",
	" for col in range(self.cols):\n",
	" for row in range(self.rows):\n",
	" result.set(row, col, round(self.get(row, col) - other.get(row, col), self.rounding))\n",
	" return result\n",
	" def multiply(self, other, *, name = None, description = None):\n",
	" # Calculate self * other\n",
	" if isinstance(other, Matrix):\n",
	" if not self.cols == other.rows:\n",
	" raise ValueError(\"Colums of left side {} have to match rows of right side {}\".format(self.cols, other.rows))\n",
	" result = Matrix(self.rows, other.cols, name = name, description = description)\n",
	" for row in range(self.rows):\n",
	" for col in range(other.cols):\n",
	" t = 0\n",
	" for k in range(self.cols):\n",
	" t = t + self.get(row, k) * other.get(k, col)\n",
	" t = round(t, 4)\n",
	" result.set(row, col, round(t, self.rounding))\n",
	" return result\n",
	" else:\n",
	" try:\n",
	" otherFloat = float(other)\n",
	" result = Matrix(self.rows, self.cols, name = name, description = description)\n",
	" for r in range(self.rows):\n",
	" for c in range(self.cols):\n",
	" result.set(r,c, round(self.get(r,c,)*otherFloat, self.rounding))\n",
	" return result\n",
	" except ValueError:\n",
	" raise ValueError(\"Matrix has to be multiplied by matrix or scalar\")"
	]
	},
	{
	"cell_type": "markdown",
	"id": "f32d5200",
	"metadata": {},
	"source": [
	"As a helper let's also define a simple identity matrix class:"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 3,
	"id": "654a576a",
	"metadata": {},
	"outputs": [],
	"source": [
	"class IdentityMatrix(Matrix):\n",
	" def __init__(self, size, *, name = None, description = None):\n",
	" super().__init__(size, size, name = name, description = description)\n",
	" for k in range(size):\n",
	" self.set(k,k,1)"
	]
	},
	{
	"cell_type": "markdown",
	"id": "f692886a",
	"metadata": {},
	"source": [
	"## The QR decomposition\n",
	"\n",
	"To perform QR decomposition using Givens rotations let's first define a class that creates the $G_{i,j}$ matrix when we also tell it the values of $x_i$ and $x_k$ in our given column where $i$ and $k$ are the row indices"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 4,
	"id": "9834c499",
	"metadata": {},
	"outputs": [],
	"source": [
	"class GivensRotationMatrix(IdentityMatrix):\n",
	" def __init__(self, size, i, j, x1, x2, *, name = None, description = None):\n",
	" super().__init__(size, name = name, description = description)\n",
	" if (i >= size) or (i < 0) or (j >= size) or (j < 0):\n",
	" raise ValueError(\"Invalid first index i or j\")\n",
	" # Calculate s and c\n",
	" s = 1.0 / math.sqrt(x1 * x1 / (x2 * x2) + 1)\n",
	" c = s * x1 / x2\n",
	" self.set(i,i,round(c, self.rounding))\n",
	" self.set(j,i,round(-1 * s, self.rounding))\n",
	" self.set(i,j,round(s, self.rounding))\n",
	" self.set(j,j,round(c, self.rounding))"
	]
	},
	{
	"cell_type": "markdown",
	"id": "b96346c6",
	"metadata": {},
	"source": [
	"Now we know that we should apply the transformations to zero out the lower triangular area starting from the lowest row in the leftmost column proceeding upwards into the first subdiagonal, then continuing the same process from bottom to top in the second column, etc."
	]
	},
	{
	"cell_type": "code",
	"execution_count": 5,
	"id": "6b90b208",
	"metadata": {},
	"outputs": [],
	"source": [
	"def QRDecomposeWithGivens(matIn, *, traceOutput = False):\n",
	" matR = matIn.copy()\n",
	" matQ = IdentityMatrix(matIn.rows)\n",
	" if traceOutput:\n",
	" print(\"Input into QR decomposition:\")\n",
	" print(matIn)\n",
	" for col in range(matIn.cols):\n",
	" for row in range(matIn.rows-1, col, -1):\n",
	" if traceOutput:\n",
	" x1 = matR.get(row-1, col)\n",
	" x2 = matR.get(row,col)\n",
	" rotationMatrix = GivensRotationMatrix(matIn.rows, row-1, row, matR.get(row-1, col), matR.get(row,col))\n",
	" matR = rotationMatrix.multiply(matR)\n",
	" matQ = matQ.multiply(rotationMatrix.transpose())\n",
	" if traceOutput:\n",
	" print(\"--- Next step ---\")\n",
	" print(\"Rotation matrix for col {}, row {}; x{}: {}, x{}: {}\".format(col, row, row-1, x1, row, x2))\n",
	" print(rotationMatrix)\n",
	" print(\"New R\")\n",
	" print(matR)\n",
	" print(\"New Q\")\n",
	" print(matQ)\n",
	"\n",
	" return (matQ, matR)"
	]
	},
	{
	"cell_type": "markdown",
	"id": "08b1cfff",
	"metadata": {},
	"source": [
	"## A simple test\n",
	"\n",
	"To verify the decomposition works let's use a simple test using a matrix that has also been found in some textbook:"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 6,
	"id": "cda3ed6a",
	"metadata": {},
	"outputs": [
	{
	"name": "stdout",
	"output_type": "stream",
	"text": [
	"0.8147\t0.0975\t0.1576\t\n",
	"0.9058\t0.2785\t0.9706\t\n",
	"0.127\t0.5469\t0.9572\t\n",
	"0.9134\t0.9575\t0.4854\t\n",
	"0.6324\t0.9649\t0.8003\t\n",
	"\n"
	]
	}
	],
	"source": [
	"testMatrix1 = Matrix(5, 3, initList = [ [ 0.8147, 0.0975, 0.1576 ], [ 0.9058, 0.2785, 0.9706 ], [ 0.1270, 0.5469, 0.9572 ], [ 0.9134, 0.9575, 0.4854 ], [ 0.6324, 0.9649, 0.8003 ] ])\n",
	"print(testMatrix1)"
	]
	},
	{
	"cell_type": "markdown",
	"id": "01628469",
	"metadata": {},
	"source": [
	"Now apply our QR decomposition:"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 7,
	"id": "72a8e3da",
	"metadata": {},
	"outputs": [
	{
	"name": "stdout",
	"output_type": "stream",
	"text": [
	"Q:\n",
	"0.4927\t-0.4807\t-0.1779\t-0.7032\t0.0\t\n",
	"0.5477\t-0.3583\t0.5778\t0.4824\t-0.0706\t\n",
	"0.0768\t0.4754\t0.6342\t-0.4316\t0.4235\t\n",
	"0.5524\t0.3391\t-0.481\t0.2769\t0.5216\t\n",
	"0.3824\t0.5473\t-0.0312\t-0.0984\t-0.7373\t\n",
	"\n",
	"R:\n",
	"1.6536\t1.1405\t1.257\t\n",
	"0.0\t0.9662\t0.6341\t\n",
	"0.0\t0.0\t0.8815\t\n",
	"0.0\t0.0\t-0.0\t\n",
	"0.0\t0.0\t0.0\t\n",
	"\n"
	]
	}
	],
	"source": [
	"q,r = QRDecomposeWithGivens(testMatrix1, traceOutput = False)\n",
	"\n",
	"print(\"Q:\")\n",
	"print(q)\n",
	"print(\"R:\")\n",
	"print(r)"
	]
	},
	{
	"cell_type": "markdown",
	"id": "107c3d02",
	"metadata": {},
	"source": [
	"And calculate the product $Q*R$ again that should yield our input (up to rounding errors). In the end also calculate the difference to make comparison easier:"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 8,
	"id": "6e07cb68",
	"metadata": {},
	"outputs": [
	{
	"name": "stdout",
	"output_type": "stream",
	"text": [
	"0.8147\t0.0975\t0.1577\t\n",
	"0.9057\t0.2785\t0.9706\t\n",
	"0.127\t0.5469\t0.957\t\n",
	"0.9134\t0.9577\t0.4854\t\n",
	"0.6323\t0.9649\t0.8002\t\n",
	"\n"
	]
	}
	],
	"source": [
	"qr = q.multiply(r)\n",
	"print(qr)"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 9,
	"id": "760a34dd",
	"metadata": {},
	"outputs": [
	{
	"name": "stdout",
	"output_type": "stream",
	"text": [
	"0.0\t0.0\t-0.0001\t\n",
	"0.0001\t0.0\t0.0\t\n",
	"0.0\t0.0\t0.0002\t\n",
	"0.0\t-0.0002\t0.0\t\n",
	"0.0001\t0.0\t0.0001\t\n",
	"\n"
	]
	}
	],
	"source": [
	"difference = testMatrix1.sub(qr)\n",
	"print(difference)"
	]
	},
	{
	"cell_type": "markdown",
	"id": "c9c91011",
	"metadata": {},
	"source": [
	"As one can see multiplying the Q and R matrices again yields the initial matrix. To verify if we have an (up to rounding errors) orthogonal matrix let's check if $Q^T * Q = \\mathfrak{1}$"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 10,
	"id": "5dcb4570",
	"metadata": {},
	"outputs": [
	{
	"name": "stdout",
	"output_type": "stream",
	"text": [
	"1.0\t0.0\t-0.0001\t-0.0001\t0.0\t\n",
	"0.0\t1.0\t-0.0002\t0.0\t-0.0\t\n",
	"-0.0001\t-0.0002\t1.0\t-0.0\t-0.0001\t\n",
	"-0.0001\t0.0\t-0.0\t0.9998\t0.0001\t\n",
	"0.0\t-0.0\t-0.0001\t0.0001\t1.0\t\n",
	"\n"
	]
	}
	],
	"source": [
	"print(q.transpose().multiply(q))"
	]
	},
	{
	"cell_type": "markdown",
	"id": "555c1134",
	"metadata": {},
	"source": [
	"Up to rounding errors this yields acceptable results as an orthogonal matrix. Thus our decomposition seems to really work as expected."
	]
	},
	{
	"cell_type": "markdown",
	"id": "ec325a34",
	"metadata": {},
	"source": [
	"## Including trace output\n",
	"\n",
	"To make verification of implementations easier the same QR decomposition can trace all intermediate steps:"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 11,
	"id": "4591c08c",
	"metadata": {},
	"outputs": [
	{
	"name": "stdout",
	"output_type": "stream",
	"text": [
	"Input into QR decomposition:\n",
	"0.8147\t0.0975\t0.1576\t\n",
	"0.9058\t0.2785\t0.9706\t\n",
	"0.127\t0.5469\t0.9572\t\n",
	"0.9134\t0.9575\t0.4854\t\n",
	"0.6324\t0.9649\t0.8003\t\n",
	"\n",
	"--- Next step ---\n",
	"Rotation matrix for col 0, row 4; x3: 0.9134, x4: 0.6324\n",
	"1\t0.0\t0.0\t0.0\t0.0\t\n",
	"0.0\t1\t0.0\t0.0\t0.0\t\n",
	"0.0\t0.0\t1\t0.0\t0.0\t\n",
	"0.0\t0.0\t0.0\t0.82217275\t0.56923806\t\n",
	"0.0\t0.0\t0.0\t-0.56923806\t0.82217275\t\n",
	"\n",
	"New R\n",
	"0.8147\t0.0975\t0.1576\t\n",
	"0.9058\t0.2785\t0.9706\t\n",
	"0.127\t0.5469\t0.9572\t\n",
	"1.111\t1.3365\t0.8546\t\n",
	"0.0\t0.2483\t0.3817\t\n",
	"\n",
	"New Q\n",
	"1.0\t0.0\t0.0\t0.0\t0.0\t\n",
	"0.0\t1.0\t0.0\t0.0\t0.0\t\n",
	"0.0\t0.0\t1.0\t0.0\t0.0\t\n",
	"0.0\t0.0\t0.0\t0.8222\t-0.5692\t\n",
	"0.0\t0.0\t0.0\t0.5692\t0.8222\t\n",
	"\n",
	"--- Next step ---\n",
	"Rotation matrix for col 0, row 3; x2: 0.127, x3: 1.111\n",
	"1\t0.0\t0.0\t0.0\t0.0\t\n",
	"0.0\t1\t0.0\t0.0\t0.0\t\n",
	"0.0\t0.0\t0.11357181\t0.99352979\t0.0\t\n",
	"0.0\t0.0\t-0.99352979\t0.11357181\t0.0\t\n",
	"0.0\t0.0\t0.0\t0.0\t1\t\n",
	"\n",
	"New R\n",
	"0.8147\t0.0975\t0.1576\t\n",
	"0.9058\t0.2785\t0.9706\t\n",
	"1.1182\t1.39\t0.9578\t\n",
	"-0.0\t-0.3916\t-0.8539\t\n",
	"0.0\t0.2483\t0.3817\t\n",
	"\n",
	"New Q\n",
	"1.0\t0.0\t0.0\t0.0\t0.0\t\n",
	"0.0\t1.0\t0.0\t0.0\t0.0\t\n",
	"0.0\t0.0\t0.1136\t-0.9935\t0.0\t\n",
	"0.0\t0.0\t0.8169\t0.0934\t-0.5692\t\n",
	"0.0\t0.0\t0.5655\t0.0646\t0.8222\t\n",
	"\n",
	"--- Next step ---\n",
	"Rotation matrix for col 0, row 2; x1: 0.9058, x2: 1.1182\n",
	"1\t0.0\t0.0\t0.0\t0.0\t\n",
	"0.0\t0.62944608\t0.77704417\t0.0\t0.0\t\n",
	"0.0\t-0.77704417\t0.62944608\t0.0\t0.0\t\n",
	"0.0\t0.0\t0.0\t1\t0.0\t\n",
	"0.0\t0.0\t0.0\t0.0\t1\t\n",
	"\n",
	"New R\n",
	"0.8147\t0.0975\t0.1576\t\n",
	"1.439\t1.2554\t1.3552\t\n",
	"-0.0\t0.6585\t-0.1513\t\n",
	"0.0\t-0.3916\t-0.8539\t\n",
	"0.0\t0.2483\t0.3817\t\n",
	"\n",
	"New Q\n",
	"1.0\t0.0\t0.0\t0.0\t0.0\t\n",
	"0.0\t0.6294\t-0.777\t0.0\t0.0\t\n",
	"0.0\t0.0883\t0.0715\t-0.9935\t0.0\t\n",
	"0.0\t0.6348\t0.5142\t0.0934\t-0.5692\t\n",
	"0.0\t0.4394\t0.356\t0.0646\t0.8222\t\n",
	"\n",
	"--- Next step ---\n",
	"Rotation matrix for col 0, row 1; x0: 0.8147, x1: 1.439\n",
	"0.49267685\t0.87021234\t0.0\t0.0\t0.0\t\n",
	"-0.87021234\t0.49267685\t0.0\t0.0\t0.0\t\n",
	"0.0\t0.0\t1\t0.0\t0.0\t\n",
	"0.0\t0.0\t0.0\t1\t0.0\t\n",
	"0.0\t0.0\t0.0\t0.0\t1\t\n",
	"\n",
	"New R\n",
	"1.6536\t1.1405\t1.257\t\n",
	"-0.0\t0.5337\t0.5305\t\n",
	"0.0\t0.6585\t-0.1513\t\n",
	"0.0\t-0.3916\t-0.8539\t\n",
	"0.0\t0.2483\t0.3817\t\n",
	"\n",
	"New Q\n",
	"0.4927\t-0.8702\t0.0\t0.0\t0.0\t\n",
	"0.5477\t0.3101\t-0.777\t0.0\t0.0\t\n",
	"0.0768\t0.0435\t0.0715\t-0.9935\t0.0\t\n",
	"0.5524\t0.3128\t0.5142\t0.0934\t-0.5692\t\n",
	"0.3824\t0.2165\t0.356\t0.0646\t0.8222\t\n",
	"\n",
	"--- Next step ---\n",
	"Rotation matrix for col 1, row 4; x3: -0.3916, x4: 0.2483\n",
	"1\t0.0\t0.0\t0.0\t0.0\t\n",
	"0.0\t1\t0.0\t0.0\t0.0\t\n",
	"0.0\t0.0\t1\t0.0\t0.0\t\n",
	"0.0\t0.0\t0.0\t-0.8445395\t0.53549326\t\n",
	"0.0\t0.0\t0.0\t-0.53549326\t-0.8445395\t\n",
	"\n",
	"New R\n",
	"1.6536\t1.1405\t1.257\t\n",
	"0.0\t0.5337\t0.5305\t\n",
	"0.0\t0.6585\t-0.1513\t\n",
	"0.0\t0.4637\t0.9256\t\n",
	"0.0\t0.0\t0.1349\t\n",
	"\n",
	"New Q\n",
	"0.4927\t-0.8702\t0.0\t0.0\t0.0\t\n",
	"0.5477\t0.3101\t-0.777\t0.0\t0.0\t\n",
	"0.0768\t0.0435\t0.0715\t0.839\t0.532\t\n",
	"0.5524\t0.3128\t0.5142\t-0.3837\t0.4307\t\n",
	"0.3824\t0.2165\t0.356\t0.3857\t-0.729\t\n",
	"\n",
	"--- Next step ---\n",
	"Rotation matrix for col 1, row 3; x2: 0.6585, x3: 0.4637\n",
	"1\t0.0\t0.0\t0.0\t0.0\t\n",
	"0.0\t1\t0.0\t0.0\t0.0\t\n",
	"0.0\t0.0\t0.81762457\t0.57575173\t0.0\t\n",
	"0.0\t0.0\t-0.57575173\t0.81762457\t0.0\t\n",
	"0.0\t0.0\t0.0\t0.0\t1\t\n",
	"\n",
	"New R\n",
	"1.6536\t1.1405\t1.257\t\n",
	"0.0\t0.5337\t0.5305\t\n",
	"0.0\t0.8054\t0.4092\t\n",
	"0.0\t-0.0\t0.8439\t\n",
	"0.0\t0.0\t0.1349\t\n",
	"\n",
	"New Q\n",
	"0.4927\t-0.8702\t0.0\t0.0\t0.0\t\n",
	"0.5477\t0.3101\t-0.6353\t0.4474\t0.0\t\n",
	"0.0768\t0.0435\t0.5415\t0.6448\t0.532\t\n",
	"0.5524\t0.3128\t0.1995\t-0.6098\t0.4307\t\n",
	"0.3824\t0.2165\t0.5131\t0.1104\t-0.729\t\n",
	"\n",
	"--- Next step ---\n",
	"Rotation matrix for col 1, row 2; x1: 0.5337, x2: 0.8054\n",
	"1\t0.0\t0.0\t0.0\t0.0\t\n",
	"0.0\t0.55238123\t0.83359161\t0.0\t0.0\t\n",
	"0.0\t-0.83359161\t0.55238123\t0.0\t0.0\t\n",
	"0.0\t0.0\t0.0\t1\t0.0\t\n",
	"0.0\t0.0\t0.0\t0.0\t1\t\n",
	"\n",
	"New R\n",
	"1.6536\t1.1405\t1.257\t\n",
	"0.0\t0.9662\t0.6341\t\n",
	"0.0\t0.0\t-0.2162\t\n",
	"0.0\t0.0\t0.8439\t\n",
	"0.0\t0.0\t0.1349\t\n",
	"\n",
	"New Q\n",
	"0.4927\t-0.4807\t0.7254\t0.0\t0.0\t\n",
	"0.5477\t-0.3583\t-0.6094\t0.4474\t0.0\t\n",
	"0.0768\t0.4754\t0.2629\t0.6448\t0.532\t\n",
	"0.5524\t0.3391\t-0.1505\t-0.6098\t0.4307\t\n",
	"0.3824\t0.5473\t0.103\t0.1104\t-0.729\t\n",
	"\n",
	"--- Next step ---\n",
	"Rotation matrix for col 2, row 4; x3: 0.8439, x4: 0.1349\n",
	"1\t0.0\t0.0\t0.0\t0.0\t\n",
	"0.0\t1\t0.0\t0.0\t0.0\t\n",
	"0.0\t0.0\t1\t0.0\t0.0\t\n",
	"0.0\t0.0\t0.0\t0.98746326\t0.15784903\t\n",
	"0.0\t0.0\t0.0\t-0.15784903\t0.98746326\t\n",
	"\n",
	"New R\n",
	"1.6536\t1.1405\t1.257\t\n",
	"0.0\t0.9662\t0.6341\t\n",
	"0.0\t0.0\t-0.2162\t\n",
	"0.0\t0.0\t0.8546\t\n",
	"0.0\t0.0\t-0.0\t\n",
	"\n",
	"New Q\n",
	"0.4927\t-0.4807\t0.7254\t0.0\t0.0\t\n",
	"0.5477\t-0.3583\t-0.6094\t0.4418\t-0.0706\t\n",
	"0.0768\t0.4754\t0.2629\t0.7207\t0.4235\t\n",
	"0.5524\t0.3391\t-0.1505\t-0.5342\t0.5216\t\n",
	"0.3824\t0.5473\t0.103\t-0.0061\t-0.7373\t\n",
	"\n",
	"--- Next step ---\n",
	"Rotation matrix for col 2, row 3; x2: -0.2162, x3: 0.8546\n",
	"1\t0.0\t0.0\t0.0\t0.0\t\n",
	"0.0\t1\t0.0\t0.0\t0.0\t\n",
	"0.0\t0.0\t-0.24525723\t0.96945804\t0.0\t\n",
	"0.0\t0.0\t-0.96945804\t-0.24525723\t0.0\t\n",
	"0.0\t0.0\t0.0\t0.0\t1\t\n",
	"\n",
	"New R\n",
	"1.6536\t1.1405\t1.257\t\n",
	"0.0\t0.9662\t0.6341\t\n",
	"0.0\t0.0\t0.8815\t\n",
	"0.0\t0.0\t-0.0\t\n",
	"0.0\t0.0\t0.0\t\n",
	"\n",
	"New Q\n",
	"0.4927\t-0.4807\t-0.1779\t-0.7032\t0.0\t\n",
	"0.5477\t-0.3583\t0.5778\t0.4824\t-0.0706\t\n",
	"0.0768\t0.4754\t0.6342\t-0.4316\t0.4235\t\n",
	"0.5524\t0.3391\t-0.481\t0.2769\t0.5216\t\n",
	"0.3824\t0.5473\t-0.0312\t-0.0984\t-0.7373\t\n",
	"\n"
	]
	}
	],
	"source": [
	"q,r = QRDecomposeWithGivens(testMatrix1, traceOutput = True)"
	]
	}
	],
	"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.8.12"
	}
	},
	"nbformat": 4,
	"nbformat_minor": 5
	}