4
$\begingroup$

This is a follow-up of a previous question.

Let $p$ be a polynomial with floating-point coefficients.

Is there a method for finding intervals where evaluating $p$ in floating-point arithmetic always gives the correct sign?

I want to find (ideally, maximal) disjoint interval $I_1, \dots, I_m$ such that for all floating-point numbers $a \in I_k$, the sign of $f(a)$ is always correct.

This is clearly related to isolating the zeros of $f$.

A standard example, where the sign of $f$ is wrong around a zero is $f(x)=x^6 - 6 x^5 + 15 x^4 - 20 x^3 + 15 x^2 - 6 x + 1$ near the multiple zero $x=1$. Note the $f(x)$ is the expansion of $(x-1)^6$.

$\endgroup$
4
$\begingroup$

Yes. You can compute a running error bound, i.e, a number $\mu$ such that the difference between the exact value of $y = p(x)$ and the computed value satisfies $\hat{y}$ satisfies $$|y - \hat{y}| \leq \mu u.$$ Here $u$ is the unit roundoff. You can trust the sign computed sign of $y$, when $|\hat{y}| > \mu u$.

Let $p(x) = \sum_{j=0}^n a_j x^j$, then Horner's method computes $$p_0 = a_n, \quad p_i = x p_{i-1} + a_{n-i}, \quad i = 1,2,\dotsc,n.$$ If $a_i$ and $x$ are machine numbers, then a running error bound can be computed as follows $$\mu_0 = 0, \quad z_j = p_{j-1} x, \quad p_j = z_j + a_{n-j}, \quad \mu_j = \mu_{j-1} |x| + |z_j| + |p_j|, \quad j=1,2,\dotsc, n.$$ The algorithm returns $y = p_n$ and $\mu = \mu_n$.

It is possible to reduce the cost of this algorithm from 5n flops to 4n flops, see Higham's book "Accuracy and Stability of Numerical Algorithms" for details.

$\endgroup$
  • $\begingroup$ Thanks! How can I use this to find the intervals I seek? $\endgroup$ – lhf Nov 5 '17 at 15:00
  • $\begingroup$ @lfh If I had to I would try to solve the equation $|y| = \mu u$ but I suspect this is harder than one might think. In the past I have only needed to know if I could trust the computed sign of specific $y = p(x)$. This is needed when searching robustly for a root. Then I need a bracket around the root. Even with a non-robust method like Newton's method there is no point in continuing the search once $|\hat{y}| \leq \mu u$, because the true value of $y$ could be zero for all that we know. This doesn't change that your question is valid and I don't have a good answer for you. $\endgroup$ – Carl Christian Nov 5 '17 at 16:37
  • $\begingroup$ Indeed, knowing when to trust the computed sign while searching robustly for a root is one main application and motivation. Thanks. $\endgroup$ – lhf Nov 5 '17 at 17:40
2
$\begingroup$

I want to add that in addition to Carl Christian's suggestion of using a running error bound, you can also take the general relative error bound $$ \frac{|\hat p(x)-p(x)|}{|p(x)|} \leq \gamma_{2n}\,\mathrm{cond}(p,x),\qquad \gamma_{2n} \approx 2nu,\\ \mathrm{cond}(p,x) = \frac{\sum |a_i||x|^i}{|\sum a_i x^i|} = \frac{\tilde p(x)}{|p(x)|},$$ (see http://www-pequan.lip6.fr/~jmc/polycopies/Compensation-horner.pdf), and then the condition that the sign is right can be imposed be imposing $$ \gamma_{2n}\mathrm{cond}(p,x) \leq 1, \quad\text{or}\quad |p(x)|\gtrsim 2nu\tilde p(x)$$

Essentially all this is doing is identifying those $x$ for which evaluating the polynomial at $x$ is well-conditioned, so it's just a bound on the condition number.

For $(x-1)^6$, this would result in something like $|x-1| > c u^{1/6}$ where $c$ is a small constant $c$ and $u$ is the unit roundoff. For a simple root $\alpha$, it would be $|x-\alpha| > c u$.

$\endgroup$

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.