{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have a square matrix $R$.  We consider the error for $T = R'R$ where $R'$ is the transpose of $R$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The elements of $R$ are $r_{i,j}$, where $i = 1 \\dots N, j = 1 \\dots N$.\n",
    "\n",
    "$r_{i, *}$ is row $i$ of $R$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let $R$ be a rotation matrix.  $T$ at infinite precision will be the identity matrix $I$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assume the maximum error in the specification of values $r_{i, j}$ is constant, $\\delta$. That is, any floating point value $r_{i, j}$ represents an infinite precision value between $r_{i, j} \\pm \\delta$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[r_0_0, r_0_1, r_0_2]\n",
       "[r_1_0, r_1_1, r_1_2]\n",
       "[r_2_0, r_2_1, r_2_2]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import Symbol, symarray, Matrix, matrices, simplify, nsimplify\n",
    "R = Matrix(symarray('r', (3,3)))\n",
    "R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[         r_0_0**2 + r_1_0**2 + r_2_0**2, r_0_0*r_0_1 + r_1_0*r_1_1 + r_2_0*r_2_1, r_0_0*r_0_2 + r_1_0*r_1_2 + r_2_0*r_2_2]\n",
       "[r_0_0*r_0_1 + r_1_0*r_1_1 + r_2_0*r_2_1,          r_0_1**2 + r_1_1**2 + r_2_1**2, r_0_1*r_0_2 + r_1_1*r_1_2 + r_2_1*r_2_2]\n",
       "[r_0_0*r_0_2 + r_1_0*r_1_2 + r_2_0*r_2_2, r_0_1*r_0_2 + r_1_1*r_1_2 + r_2_1*r_2_2,          r_0_2**2 + r_1_2**2 + r_2_2**2]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T = R.T * R\n",
    "T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the same result with error $\\delta$ added to each element"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[d + r_0_0, d + r_0_1, d + r_0_2]\n",
       "[d + r_1_0, d + r_1_1, d + r_1_2]\n",
       "[d + r_2_0, d + r_2_1, d + r_2_2]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = Symbol('d')\n",
    "E = matrices.ones((3,3)) * d\n",
    "RE = R + E\n",
    "RE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the result $T$ with error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[                           (d + r_0_0)**2 + (d + r_1_0)**2 + (d + r_2_0)**2, (d + r_0_0)*(d + r_0_1) + (d + r_1_0)*(d + r_1_1) + (d + r_2_0)*(d + r_2_1), (d + r_0_0)*(d + r_0_2) + (d + r_1_0)*(d + r_1_2) + (d + r_2_0)*(d + r_2_2)]\n",
       "[(d + r_0_0)*(d + r_0_1) + (d + r_1_0)*(d + r_1_1) + (d + r_2_0)*(d + r_2_1),                            (d + r_0_1)**2 + (d + r_1_1)**2 + (d + r_2_1)**2, (d + r_0_1)*(d + r_0_2) + (d + r_1_1)*(d + r_1_2) + (d + r_2_1)*(d + r_2_2)]\n",
       "[(d + r_0_0)*(d + r_0_2) + (d + r_1_0)*(d + r_1_2) + (d + r_2_0)*(d + r_2_2), (d + r_0_1)*(d + r_0_2) + (d + r_1_1)*(d + r_1_2) + (d + r_2_1)*(d + r_2_2),                            (d + r_0_2)**2 + (d + r_1_2)**2 + (d + r_2_2)**2]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TE = RE.T * RE\n",
    "TE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Subtract the true result to get the absolute error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[                  d*(3*d + 2*r_0_0 + 2*r_1_0 + 2*r_2_0), d*(3*d + r_0_0 + r_0_1 + r_1_0 + r_1_1 + r_2_0 + r_2_1), d*(3*d + r_0_0 + r_0_2 + r_1_0 + r_1_2 + r_2_0 + r_2_2)]\n",
       "[d*(3*d + r_0_0 + r_0_1 + r_1_0 + r_1_1 + r_2_0 + r_2_1),                   d*(3*d + 2*r_0_1 + 2*r_1_1 + 2*r_2_1), d*(3*d + r_0_1 + r_0_2 + r_1_1 + r_1_2 + r_2_1 + r_2_2)]\n",
       "[d*(3*d + r_0_0 + r_0_2 + r_1_0 + r_1_2 + r_2_0 + r_2_2), d*(3*d + r_0_1 + r_0_2 + r_1_1 + r_1_2 + r_2_1 + r_2_2),                   d*(3*d + 2*r_0_2 + 2*r_1_2 + 2*r_2_2)]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TTE = TE-T\n",
    "TTE.simplify()\n",
    "TTE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(d*(3*d + 2*r_0_0 + 2*r_1_0 + 2*r_2_0),\n",
       " d*(3*d + 2*r_0_1 + 2*r_1_1 + 2*r_2_1),\n",
       " d*(3*d + 2*r_0_2 + 2*r_1_2 + 2*r_2_2))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TTE[0,0], TTE[1,1], TTE[2,2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming $\\delta$ is small ($\\delta^2$ is near zero) then the diagonal values $TTE_{k, k}$ are approximately $2\\delta \\Sigma_i r_{i, k}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\Sigma_i r_{i, k}$ is the column sum for column $k$ of $R$.  We know that the column $L^2$ norms of $R$ are each 1.  We know $\\|x\\|_1 \\leq \\sqrt{n}\\|x\\|_2$ - https://en.wikipedia.org/wiki/Lp_space. Therefore the column sums must be $\\le \\sqrt{N}$.  Therefore the maximum error for the diagonal of $T$ is $\\sqrt{N} 2\\delta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More generally, the elements $k, m$ of $TTE$ are approximately $\\delta (\\Sigma_i{r_{i, k}} + \\Sigma_i{r_{i, m}})$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the error for each of the elements of $TTE$ is also bounded by $\\sqrt{N} 2\\delta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now consider the floating point calculation error.  This depends on the floating point representation we use for the calculations. Let $\\epsilon = x-1$ where $x$ is the smallest number greater than 1 that is representable in our floating point format (see https://matthew-brett.github.io/pydagogue/floating_error.html). The largest error for a calculation resulting in a value near 1 is $\\frac{\\epsilon}{2}$. For the diagonal values, the calculation error will be the error for the $r_{i,*} r_{i, *}'$ dot product.  This comprises $N$ scalar products with results each bounded by 1 ($r_{i, j} r_{i, j}$) followed by $N-1$ sums each bounded by 1.  Maximum error is therefore $(2N-1) \\frac{\\epsilon}{2}$ = $\\frac{5}{2} \\epsilon$ where $N=3$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the off-diagonal values $T_{k, m}$, we have the $r_{k,*} r_{m, *}'$ dot product.  Because $R$ is a rotation matrix, by definition this result must be zero.\n",
    "\n",
    "Because the column and row $L^2$ norms of $R$ are 1, the values in $R$ cannot be greater than 1.  Therefore $r_{k,*} r_{m, *}'$ consists of the $N$ products with results each bounded by 1 ($r_{k, j} r_{m, j}$) followed by $N-2$ sums each bounded by 1 (the last of the sums must be approximately 0).  Maximum error is therefore $(2N-2) \\frac{\\epsilon}{2}$ = $2\\epsilon$ where $N=3$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, assuming an initial error of $\\delta$ per element, and $N=3$, the maximum error for diagonal elements is $\\sqrt{3} 2 \\delta + \\frac{5}{3} \\epsilon$. For the off-diagonal elements it is $\\sqrt{3} 2 \\delta + 2 \\epsilon$."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
