6
$\begingroup$

I have two sets of large positive numbers $a_1,\ldots,a_n$ and $b_1,\ldots,b_n$. By 'large' I mean of the order of $10^{10}$. I want to calculate the ratio $$R = \frac{a_1 - a_2 + \cdots +(-1)^{n+1}a_n}{b_1 - b_2 + \cdots +(-1)^{n+1}b_n}.$$ It is already known that the exact value of $R$ is between $0$ and $1$. Each of the numerator and denominator is expected to be of the order of $10^{-1}$ to $10^1$.

I would like to know if there is an algorithm to compute ratios like $R$ stably (exactly or approximately), so that the numbers involved in each step are not too large to cause overflow errors using finite-precision arithmetic. I am using the usual 64-bit double data type in Matlab.

I am aware of the Kahan summation algorithm, but using it to compute the numerator and denominator separately may not guarantee that their ratio is also as accurate.

$\endgroup$
  • 1
    $\begingroup$ There are multiple ways in which division itself is done numerically and some of them employ subtraction as the fundamental operation - restorative division being one of them. You can use Kahan summation on even the division algorithm and maintain the same level of accuracy. More information on web.stanford.edu/class/ee486/doc/chap5.pdf $\endgroup$ – gpavanb Aug 1 '17 at 21:25
  • 3
    $\begingroup$ Why is the division the problem here? Doesn't the alternation in the enumerator and denominator already take care of the problem and lead to moderately sized terms for the division? I would have imagined that any round-off problem appears in the sums themselves. $\endgroup$ – Wolfgang Bangerth Aug 6 '17 at 22:02
1
$\begingroup$

Using Kahan summation algorithm should suffice your purpose. Having: $$ R=\frac{\sum\limits_{i=1}^{N}a_i}{\sum\limits_{i=1}^{N}b_i}=\frac{A}{B} $$ if $A$ and $B$ are computed accurate enough, the final ratio $R$ should not be a problem.

If $A$ and $B$ are computed using Kahan summation, you can expect the relative error to be $\mathcal{O}(\epsilon_\text{mach}+N\epsilon_\text{mach}^2)$. For double-precision and $N\approx 10^{10}$ you are safe ($N\ll 1/\epsilon_\text{mach}$), and your error is practically independent of $N$: $\mathcal{O}(\epsilon_\text{mach})$.

It's worth to mention, that the constant in front of the error estimate (condition # of summation $\kappa_{\sum}$) can still be relatively large. If summation condition number:

$$ \kappa_{\sum_A}=\frac{\sum\limits_{i=1}^{N}|a_i|}{\left|\sum\limits_{i=1}^{N}a_i\right|} $$

is large, your best bets would probably be in Shewchuk's algorithm. With that, you will be able to calculate your numerator and denominator to full double precision.

If $\tilde{A}$ and $\tilde{B}$ are calculated in full double precision (say using Shewchuk's algorithm): $$ \tilde{A}=A(1+\delta_A),\quad \tilde{B}=B(1+\delta_B),\quad |\delta_A|,|\delta_B|\le\epsilon_\text{mach} $$ , then

$$ \begin{aligned} \tilde{R}&=\tilde{A}\oslash\tilde{B}&=\frac{A(1+\delta_A)}{B(1+\delta_B)}(1+\delta_{\oslash})=\frac{A}{B}(1+\theta)\\ &&|\theta|\leq\frac{3\epsilon_\text{mach}}{1-3\epsilon_\text{mach}} \end{aligned} $$ Here, $\tilde{A},\tilde{B},\tilde{R}$ are floating point (IEEE-754) representation of $A$, $B$, and $R$, respectively, $\oslash$ - floating point analog of division.

Now, even if $A$ and $B$ are not calculated to the full precision (Kahan), it will just increase $\delta_A$ and $\delta_B$ (not allowing approximation by $\theta$ with such a tight bound). However, I still do not see how the division can be the problem here. Especially, since the true expected values of $A$ and $B$ are in $[10^{-1}, 10^{1}]$.

$\endgroup$
  • 2
    $\begingroup$ I think there is some confusion on the word "large". I think the OP means that his numbers $a_i$ and $b_i$ are large, not that $N$ is $\mathcal{O}(10^{10})$. So I don't know why he just cannot scale his numbers (if they are all the same order of magnitude). $\endgroup$ – GertVdE Apr 20 '18 at 8:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.