Computing Validated Running Error Bounds

The problem with the previous error bounds is that they are in terms of quantities like $\sum|a_{i}|$, which are not known in advance, and hard to compute (well, you certainly wouldn’t want to compute them with floating point arithmetic!).

For Horner’s Rule, there is an algorithm that helps:

Let $r=a_{n}$
Let $s=o(|a_{n}|/2)$
For i = n-1 down to 1 (inclusive?):
- $r=o(o(rx)+a_{i})$
- $s=o(o(s|x|)+|r|)
- (This is Horner’s rule, more or less)
$r=o(o(rx)+a_{0})$
$b=o(2o(s|x|)+|r|)$
$b=\mathbf{u}o(b/(1-(2n-1)\mathbf{u}))$

Now $r$ contains the value of the polynomial at $x$ (with errors), and $b$ is an upper bound on the error. Both are floating point numbers.

Note: One can use FMA in the algorithm above.

The number of flops for this is $4n+1$ (vs $2n$ for the plain Horner’s Rule). The reason for the more flops is the addition of absolute value and factors of 2.

Now for proof that the algorithm works.

Let $r_{i}$ be the value of $r$ after the $i$ iteration. Then $r_{i}=o(o(r_{i+1}x)+a_{i})$. (Remember we are going backwards) This is for $i<n$. For $n,r_{n}=a_{n}$. Using the standard models, we have:

\begin{equation*} (1+\epsilon_{i})r_{i}=r_{i+1}x(1+\delta_{i})+a_{i},|\epsilon_{i}|,|\delta_{i}|<\mathbf{u} \end{equation*}

(Note I used both of the equations in the standard model to get this).

Define $q_{i}=\sum_{h=i}^{n}a_{h}x^{h-i}$. This is the “correct” value of $r_{i}$. Define $e_{i}=r_{i}-q_{i}$. Note that $e_{n}=0$. Also note that $q_{i}=q_{i+1}x+a_{i}$. Using these, we can show:

\begin{equation*} e_{i}=xe_{i+1}+\delta_{i}xr_{i+1}-\epsilon_{i}r_{i} \end{equation*}

Taking absolute values:

\begin{equation*} e_{i}\le |x||e_{i+1}|+\mathbf{u}(|x||r_{i+1}|+|r_{i}|) \end{equation*}

Now:

\begin{equation*} \left|r-\sum_{i=0}^{n}a_{i}x^{i}\right|=|e_{0}|\le \mathbf{u}E_{0} \end{equation*}

$E_{0}=0$ and $E_{i}=(E_{i+1}+|r_{i+1}|)|x|+|r_{i}|$. This leads to:

\begin{equation*} E_{0}=|r_{n}x^{n}|+2\sum_{i=1}^{n-1}|r_{i}x^{i}|+|r_{0}| \end{equation*}

\begin{equation*} E_{0}=2S(|x|)|x|+|r| \end{equation*}

where

\begin{equation*} S(x)=\frac{|a_{n}|}{2}x^{n-1}+\sum_{i=0}^{n-2}|r_{i+1}|x^{i} \end{equation*}

This is a polynomial of degree $n-1$, and has non-negative coefficients. Therefore, $S(|x|)\le s(1+\mathbf{u})^{2n-2}$ where $s$ is the floating point number you get from the algorithm. Then use the last theorem in the previous section to get the conclusion.

To be frank, I got bored halfway and didn’t verify the last few steps.