We have a square matrix $R$.  We consider the error for $T = R'R$ where $R'$ is the transpose of $R$.

The elements of $R$ are $r_{i,j}$, where $i = 1 \dots N, j = 1 \dots N$.

$r_{i, *}$ is row $i$ of $R$.

Now let $R$ be a rotation matrix.  $T$ at infinite precision will be the identity matrix $I$

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$


```
from sympy import Symbol, symarray, Matrix, matrices, simplify, nsimplify
R = Matrix(symarray('r', (3,3)))
R
```




    [r_0_0, r_0_1, r_0_2]
    [r_1_0, r_1_1, r_1_2]
    [r_2_0, r_2_1, r_2_2]




```
T = R.T * R
T
```




    [         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]
    [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]
    [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]



Now the same result with error $\delta$ added to each element


```
d = Symbol('d')
E = matrices.ones((3,3)) * d
RE = R + E
RE
```




    [d + r_0_0, d + r_0_1, d + r_0_2]
    [d + r_1_0, d + r_1_1, d + r_1_2]
    [d + r_2_0, d + r_2_1, d + r_2_2]



Calculate the result $T$ with error


```
TE = RE.T * RE
TE
```




    [                           (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)]
    [(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)]
    [(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]



Subtract the true result to get the absolute error


```
TTE = TE-T
TTE.simplify()
TTE
```




    [                  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)]
    [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)]
    [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)]




```
TTE[0,0], TTE[1,1], TTE[2,2]
```




    (d*(3*d + 2*r_0_0 + 2*r_1_0 + 2*r_2_0),
     d*(3*d + 2*r_0_1 + 2*r_1_1 + 2*r_2_1),
     d*(3*d + 2*r_0_2 + 2*r_1_2 + 2*r_2_2))



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}$

$\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$ - http://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$.

More generally, the elements $k, m$ of $TTE$ are approximately $\delta (\Sigma_i{r_{i, k}} + \Sigma_i{r_{i, m}})$

So the error for each of the elements of $TTE$ is also bounded by $\sqrt{N} 2\delta$.

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 http://matthew-brett.github.com/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$.

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.

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$.

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$.

cf35cda2bee9298d5a963b9e546723741eeb8ee2
ipynb
Jupyter Notebook
doc/source/notebooks/ata_error.ipynb
tobon/nibabel
ff2b5457207bb5fd6097b08f7f11123dc660fda7
BSD-3-Clause
2015-10-01T01:13:59.000Z
2015-10-01T01:13:59.000Z
doc/source/notebooks/ata_error.ipynb
tobon/nibabel
ff2b5457207bb5fd6097b08f7f11123dc660fda7
BSD-3-Clause
2015-11-13T03:05:24.000Z
2016-08-06T19:18:54.000Z
doc/source/notebooks/ata_error.ipynb
tobon/nibabel
ff2b5457207bb5fd6097b08f7f11123dc660fda7
BSD-3-Clause
2019-02-27T20:48:03.000Z
2019-02-27T20:48:03.000Z
Qwen/Qwen-72B
1. YES
2. YES