Backward Error Measures for Roots of Polynomials

In this article we analyze the problem of measuring the backward error for a set of approximations for the roots of a polynomial with complex coefficients. For a general introduction to the notion of backward error analysis, the reader can consult for instance [Hig02, Section 1.5]. Consider a set of approximate solutions $\hat{X}=\{\hat{x}_{1},\ldots,\hat{x}_{d}\}\subset\mathbb{C}^{*}=\mathbb{C}\setminus\{0\}$ of a polynomial equation with nonzero coefficients

$$f=c_{0}+c_{1}x+\ldots+c_{d-1}x^{d-1}+c_{d}x^{d}=c_{d}(x-x_{1})\cdots(x-x_{d})=0,\quad f\in\mathbb{C}[x],$$

with solutions $X=\{x_{1},\ldots,x_{d}\}\subset\mathbb{C}^{*}$. The backward error of $\hat{X}$ is a measure for the ‘distance’ of $f$ to the polynomial

$$\hat{f}=\hat{c}_{0}+\hat{c}_{1}x+\ldots+\hat{c}_{d-1}x^{d-1}+c_{d}x^{d}=c_{d}(x-\hat{x}_{1})\cdots(x-\hat{x}_{d}),$$

whose roots are exactly the points in $\hat{X}$. How to measure this distance turns out to be a surprisingly subtle problem. A first and natural measure is the 2-norm distance between the coefficients of $f$ and $\hat{f}$. The norm-wise backward error (NBE) of $\hat{X}$ is

$$\textup{NBE}(\hat{X})=\sqrt{\frac{|c_{0}-\hat{c}_{0}|^{2}+\cdots+|c_{d-1}-\hat{c}_{d-1}|^{2}}{|c_{0}|^{2}+\cdots+|c_{d}|^{2}}}=\frac{\lVert{c-\hat{c}}\rVert_{2}}{\lVert{c}\rVert_{2}},$$

where $c=(c_{0},\ldots,c_{d})\in\mathbb{C}^{d},\hat{c}=(\hat{c}_{0},\ldots,c_{d})\in\mathbb{C}^{d}$. In [AMVW15], an algorithm is proposed that computes a set of approximate solutions $\hat{X}$ satisfying $\textup{NBE}(\hat{X})=O(u)$, where $u$ is the unit round-off. Such an algorithm is called norm-wise backward stable. However, it turns out that this type of stability is too ‘weak’ in a sense we explain by means of an example. Consider the polynomial $f=a(x-10^{6})(x-10^{-6})$. The set of approximate solutions $\hat{X}=\{10^{6}+u10^{6},10^{-6}+u\}$ would satisfy $\textup{NBE}(\hat{X})=O(u)$. Indeed,

$$\hat{c}=a(1+(u10^{6}+u+u^{2}10^{6}),-(10^{6}+10^{-6})-u(10^{6}+1),1),\quad c=a(1,-(10^{6}+10^{-6}),1)$$

and $\lVert{c-\hat{c}}\rVert_{2}/\lVert{c}\rVert_{2}=O(u)$. This means that we would allow a relative error of size $u10^{6}$ on the coefficient vector $c$. However, computing the roots of $f$ with $a=0.2$ using the Julia package PolynomialRoots (using the command roots) we get

julia> abs.((c - chat)./c)
3-element Array{Float64,1}:
 2.7755575615628914e-16
 0.0
 0.0

which shows that we can obtain better element-wise accuracy on the coefficient vector. This suggests another, more ‘strict’, measure for the backward error. The element-wise backward error (EBE) of $\hat{X}$ is

$$\textup{EBE}(\hat{X})=\max_{i=0,\ldots,d-1}\left|\frac{c_{i}-\hat{c}_{i}}{c_{i}}\right|.$$

Unfortunately, this measure turns out to be too strict. In [MVD15], Mastronardi and Van Dooren show that no algorithm for finding the roots of a general quadratic polynomial is element-wise backward stable, meaning that it computes $\hat{X}$ such that $\textup{EBE}(\hat{X})=O(u)$. As an alternative measure, the authors of [MVD15] propose the following definition in the case where $d=2$. The element-wise mixed backward error (EMBE) of $\hat{X}$, denoted $\textup{EMBE}(\hat{X})$, is the smallest number $\varepsilon\geq 0$ such that there exists some $\tilde{X}=\{\tilde{x}_{1},\ldots,\tilde{x}_{d}\}\subset\mathbb{C}$ and

$$\tilde{f}=\tilde{c}_{0}+\tilde{c}_{1}x+\cdots+\tilde{c}_{d-1}x^{d-1}+c_{d}x^{d}=c_{d}(x-\tilde{x}_{1})\cdots(x-\tilde{x}_{d})$$

such that

	$\displaystyle|\hat{x}_{i}-\tilde{x}_{i}|$	$\displaystyle\leq\varepsilon|\tilde{x}_{i}|,\quad i=1,\ldots,d,$
	$\displaystyle|c_{i}-\tilde{c}_{i}|$	$\displaystyle\leq\varepsilon|c_{i}|,\quad i=0,\ldots,d-1.$

Note that the second set of inequalities is equivalent to $\textup{EBE}(\tilde{X})\leq\varepsilon$. In the same paper, the authors also show the implication $\textup{EMBE}(\hat{X})=O(u)\Rightarrow\textup{NBE}(\hat{X})=O(u)$ in the case where $d=2$. The implication $\textup{EBE}(\hat{X})=O(u)\Rightarrow\textup{EMBE}(\hat{X})=O(u)$ is obvious from the definition.
The advantage of EMBE as a measure for the backward error is that it results in a notion of stability, i.e. element-wise mixed backward stability, which is stronger than norm-wise backward stability and can provably be obtained for quadratic polynomials [MVD15]. A drawback of this measure is that it is hard to compute $\textup{EMBE}(\hat{X})$ for a given set of approximate solutions $\hat{X}$ because of the rather abstract definition. In [TVB20], Tisseur and Van Barel define the min-max element-wise backward error of $\hat{X}$ as

$$\textup{TBE}(\hat{X})=\max_{i=0,\ldots,d-1}\left|\frac{c_{i}-\hat{c}_{i}}{r_{i}c_{i}}\right|$$

where $r_{i}\geq 1$ are constants that can be computed in linear time from the coefficients of $f$. The $r_{i}$ depend only on the tropical roots of $f$, which is why we will refer to this error measure as the tropical backward error (TBE). We will give a definition of the numbers $r_{i}$ in Section 3. The authors of [TVB20] also provide an algorithm that, under some assumptions on the numerical behavior of a modified QZ-algorithm (see [TVB20, Section 5, Assumption 1]), computes a set of approximate roots $\hat{X}$ satisfying $\textup{TBE}(\hat{X})=O(u)$.

In this paper, we investigate the relations between the TBE and the EMBE. In particular, we show that these error measures are equivalent under suitable assumptions. Here’s a simplified version of our first main theorem.

The strategy for proving Theorem 1.1 will also allow us to prove that, under the assumptions of the theorem, $\textup{EMBE}(\hat{X})=O(u)$ implies $\textup{NBE}(\hat{X})=O(u)$. This was proved in [MVD15] for the case where $d=2$. Under some stronger assumptions, we also prove the reverse implication.

We will give numerical evidence that Theorem 1.2 holds without the extra assumption on the polynomial $f$. In summary, we have the following diagram, where the arrows are implications.

(1)

Here, the black arrows are implications which are obvious from the definitions. The blue arrows are implications that we prove under the assumptions of Theorem 1.1. The dashed arrow represents Theorem 1.2, which uses stronger assumptions.

This article is organized as follows. In the next section we prove the equivalence of the tropical and element-wise mixed backward error measures in the case where $d=2$. This can be seen as an extension of the analysis in [MVD15], and it is instructive for the general case. The proofs for general $d$ are given in Section 3. In Section 4 we show some computational experiments and give numerical evidence that Theorem 1.2 holds under weaker assumptions.

In this section, we prove that the element-wise mixed backward error (EMBE) as introduced by Mastronardi and Van Dooren [MVD15] and the tropical backward error (TBE) from [TVB20] are equivalent backward error measures for the roots of a quadratic polynomial

$$f=ax^{2}+bx+c=a(x-x_{1})(x-x_{2}).$$

For simplicity, we assume that $a,b,c\in\mathbb{C}^{*}=\mathbb{C}\setminus\{0\}$. For the approximate roots $\hat{X}=\{\hat{x}_{1},\hat{x}_{2}\}$ of $f$, $\textup{EMBE}(\hat{X})$ is the smallest number $\varepsilon\geq 0$ such that there exists $\tilde{X}=\{\tilde{x}_{1},\tilde{x}_{2}\}$ with

	$\displaystyle\tilde{f}=a(x-\tilde{x}_{1})(x-\tilde{x}_{2})$	$\displaystyle=ax^{2}+\tilde{b}x+\tilde{c},$
	$\displaystyle|\hat{x}_{1}-\tilde{x}_{1}|\leq\varepsilon|\tilde{x}_{1}|,\quad$	$\displaystyle|\hat{x}_{2}-\tilde{x}_{2}|\leq\varepsilon|\tilde{x}_{2}|,$
	$\displaystyle|b-\tilde{b}|\leq\varepsilon|b|,\quad$	$\displaystyle|c-\tilde{c}|\leq\varepsilon|c|.$

In analogy with [TVB20], we define the $\textup{TBE}(\hat{X})$ to be the smallest number $\varepsilon\geq 0$ such that

	$\displaystyle\hat{f}=a(x-\hat{x}_{1})(x-\hat{x}_{2})=ax^{2}+\hat{b}x+\hat{c}$
	$\displaystyle|b-\hat{b}|\leq r_{b}\varepsilon|b|,\quad|c-\hat{c}|\leq\varepsilon|c|,$

where $r_{b}=\max(1,\sqrt{|ac|}/|b|)$. The definition of $r_{b}$ will be clarified in Section 3. For the approximate roots $\hat{x}_{j}$ and $\tilde{x}_{j}$ in these definitions, we will assume that the order of magnitude of $|\tilde{x}_{j}|$ and $|\hat{x}_{j}|$ is the same as the order of magnitude of $|x_{j}|$ (that is, we allow relative errors of size at most 1). Note that by the definition of EMBE, it is sufficient that this is satisfied for $|\hat{x}_{j}|$.

We will now relate these two error measures. Let $\varepsilon=\textup{EMBE}(\hat{X})$. We observe

	$\displaystyle\frac{\hat{c}-c}{c}$	$\displaystyle=\frac{\hat{c}-\tilde{c}}{c}+\frac{\tilde{c}-c}{c}=\frac{\hat{x}_{1}\hat{x}_{2}-\tilde{x}_{1}\tilde{x}_{2}}{x_{1}x_{2}}+\frac{\tilde{c}-c}{c}$
		$\displaystyle=\frac{(\tilde{x}_{1}+(\hat{x}_{1}-\tilde{x}_{1}))(\tilde{x}_{2}+(\hat{x}_{2}-\tilde{x}_{2}))-\tilde{x}_{1}\tilde{x}_{2}}{x_{1}x_{2}}+\frac{\tilde{c}-c}{c}\,.$

Since $|c-\tilde{c}|\leq\varepsilon|c|$, we have $|1-\frac{\tilde{x}_{1}\tilde{x}_{2}}{x_{1}x_{2}}|\leq\varepsilon$ and it follows

$$\left|\frac{\hat{c}-c}{c}\right|\leq(2\varepsilon+\varepsilon^{2})\left|\frac{\tilde{x}_{1}\tilde{x}_{2}}{x_{1}x_{2}}\right|+\varepsilon\lesssim 3\varepsilon\,.$$

For the coefficient $\hat{b}$, we find in an analogous way that

$$\left|\frac{\hat{b}-b}{b}\right|\leq\varepsilon\left(1+\frac{|\tilde{x}_{1}|+|\tilde{x}_{2}|}{|x_{1}+x_{2}|}\right).$$

(2)

If the solutions $x_{1}$ and $x_{2}$ have different orders of magnitude, there does not occur any cancellation in the denominator of the right hand side of (2) and this implies $\left|\frac{\hat{b}-b}{b}\right|\lesssim 2\varepsilon$. However, if the order of magnitude of both solutions is the same, the factor standing with $\varepsilon$ may be significantly larger than 1 due to cancellation. We will now make this precise and relate this to the number $r_{b}$. We define $\gamma=x_{1}/x_{2}$. By the assumption that $|\tilde{x}_{i}|$ is of the same order of magnitude as $|x_{i}|$, (2) can be written as

$$\left|\frac{\hat{b}-b}{b}\right|\leq\varepsilon\left(1+K\frac{|\gamma|+1}{|\gamma+1|}\right)$$

(3)

with $K$ a small constant. We assume, without loss of generality, that $0<|\gamma|\leq 1$. Note that we have for $0<|\gamma|<1$ the inequality

$$\frac{|\gamma|+1}{|\gamma+1|}\leq\frac{|\gamma|+1}{||\gamma|-1|}=\frac{|\gamma|+1}{1-|\gamma|},$$

(4)

which follows from $|x-y|\geq||x|-|y||,\forall x,y\in\mathbb{C}$ applied to $x=\gamma,y=-1$. Assume that

$$\frac{\sqrt{|ac|}}{|b|}=\frac{\sqrt{|\gamma|}}{|\gamma+1|}\leq 1.$$

In this case, we have

$$\frac{|\gamma|+1}{|\gamma+1|}\leq\frac{|\gamma|+1}{\sqrt{|\gamma|}}.$$

(5)

Now, assume

$$\frac{\sqrt{|ac|}}{|b|}=\frac{\sqrt{|\gamma|}}{|\gamma+1|}\geq 1.$$

In this case

$$\frac{|\gamma|+1}{|\gamma+1|}\left(\frac{\sqrt{|ac|}}{|b|}\right)^{-1}\leq\frac{|\gamma|+1}{|\gamma+1|}\leq\frac{|\gamma|+1}{1-|\gamma|}.$$

(6)

Also,

$$\frac{|\gamma|+1}{|\gamma+1|}\left(\frac{\sqrt{|ac|}}{|b|}\right)^{-1}=\frac{|\gamma|+1}{|\gamma+1|}\left(\frac{\sqrt{|\gamma|}}{|\gamma+1|}\right)^{-1}=\frac{|\gamma|+1}{\sqrt{|\gamma|}}.$$

(7)

Using (4)-(7), we find that

$$\frac{|\gamma|+1}{|\gamma+1|}r_{b}^{-1}\leq\min\left(\frac{|\gamma|+1}{1-|\gamma|},\frac{|\gamma|+1}{\sqrt{|\gamma|}}\right),$$

which gives

$$\frac{|\gamma|+1}{|\gamma+1|}r_{b}^{-1}\leq\begin{cases}\frac{\alpha+1}{1-\alpha}=\sqrt{5}&0<|\gamma|\leq\alpha\\ \frac{\alpha+1}{\sqrt{\alpha}}=\sqrt{5}&\alpha\leq|\gamma|\leq 1\end{cases}\quad\Rightarrow\leavevmode\nobreak\ \frac{|\gamma|+1}{|\gamma+1|}r_{b}^{-1}\leq\sqrt{5},$$

where $\alpha=\frac{3}{2}-\frac{\sqrt{5}}{2}$. This is illustrated in Figure 1.

It follows immediately from this observation and (3) that

$$\left|\frac{\hat{b}-b}{b}\right|\leq\varepsilon\left(\sqrt{5}K+r_{b}^{-1}\right)r_{b}\leq\varepsilon(\sqrt{5}K+1)r_{b}.$$

Figure 2 illustrates the values of $r_{b}$ and $\frac{|\gamma|+1}{|\gamma+1|}r_{b}^{-1}$ as a function of $\gamma$ in the unit disk.

This shows that element-wise mixed backward stability implies tropical backward stability. The converse also holds, as we will now show. Let $\textup{TBE}(\hat{X})=\varepsilon$ for the computed roots $\hat{X}=\{\hat{x}_{1},\hat{x}_{2}\}$. If $r_{b}=1$ then we can take $\tilde{x}_{1}=\hat{x}_{1}$ and $\tilde{x}_{2}=\hat{x}_{2}$ such that $\textup{EMBE}(\hat{X})\leq\textup{TBE}(\hat{X})$. Suppose now that $r_{b}=\sqrt{|ac|}/|b|>1$ and without loss of generality assume $|\hat{x}_{1}|\leq|\hat{x}_{2}|$. $\textup{TBE}(\hat{X})=\varepsilon$ implies that there exists $\hat{\delta}_{b}\in\mathbb{C}$ with $|\hat{\delta}_{b}|\leq\varepsilon$ such that

$$-\hat{b}=a(\hat{x}_{1}+\hat{x}_{2})=-b(1+r_{b}\hat{\delta}_{b})\,.$$

Let $\tilde{x}_{1}=\hat{x}_{1}+\frac{b}{a}r_{b}\hat{\delta}_{b}$ and $\tilde{x}_{2}=\hat{x}_{2}$. Then we have

$$-\tilde{b}=a(\tilde{x}_{1}+\tilde{x}_{2})=a(\hat{x}_{1}+\frac{b}{a}r_{b}\hat{\delta}_{b}+\hat{x}_{2})=a(\hat{x}_{1}+\hat{x}_{2})+br_{b}\hat{\delta}_{b}=-b(1+r_{b}\hat{\delta}_{b})+br_{b}\hat{\delta}_{b}=-b\,.$$

Note that

$$|\tilde{x}_{1}-\hat{x}_{1}|=\left|\frac{b}{a}r_{b}\hat{\delta}_{b}\right|=\sqrt{\left|\frac{c}{a}\right|}|\hat{\delta}_{b}|=\sqrt{|x_{1}x_{2}|}|\hat{\delta}_{b}|\lesssim|\tilde{x}_{1}|\varepsilon\,.$$

We also have

$$\tilde{c}-\hat{c}=a\tilde{x}_{1}\tilde{x}_{2}-a\hat{x}_{1}\hat{x}_{2}=a(\hat{x}_{1}+\frac{b}{a}r_{b}\hat{\delta}_{b})\hat{x}_{2}-a\hat{x}_{1}\hat{x}_{2}=br_{b}\hat{\delta}_{b}\hat{x}_{2}$$

from which we get

$$|\tilde{c}-\hat{c}|=|a|\left|\frac{b}{a}r_{b}\hat{\delta}_{b}\hat{x}_{2}\right|\lesssim\varepsilon|a||\tilde{x}_{1}\tilde{x}_{2}|=\varepsilon|\tilde{c}|\,.$$

We conclude that tropical backward stability implies element-wise mixed backward stability.

We now generalize the results of the previous section to polynomials of arbitrary degree $d$. In the following let $f=\sum_{i=0}^{d}c_{i}x^{i}\in\mathbb{C}[x]$, $c_{0},c_{d}\neq 0$, be a polynomial of degree $d$ with roots $X=\{x_{1},\ldots,x_{d}\}\subset\mathbb{C}^{*}$ and let $\hat{X}=\{\hat{x}_{1},\ldots,\hat{x}_{d}\}\subset\mathbb{C}^{*}$ be the approximate roots. Without loss of generality we assume that the roots are labeled such that $|x_{1}|\leq|x_{2}|\leq\ldots\leq|x_{d}|$.

We first formally generalize the notion of the element-wise mixed backward error due to Mastronardi and Van Dooren [MVD15] from quadratic polynomials to general polynomials of degree $d$.

Before we can state the generalization of the tropical backward error for degree $d$ polynomials we need to introduce some definitions and concepts. We define the tropical polynomial $\texttt{t}{f}(\tau)$ associated to $f(x)$ as

$$\texttt{t}{f}:\mathbb{R}\cup\{-\infty\}\rightarrow\mathbb{R}\cup\{-\infty\},\quad\tau\mapsto\max_{0\leq i\leq d}v_{i}+i\tau$$

(8)

where $v_{i}=\log|c_{i}|$. The number $v_{i}$ is the valuation of $c_{i}$ under the valuation map $\log|\cdot|$. Any base for the logarithm can be used in theory. We want to think of the image under the valuation map as the ‘order of magnitude’ of the modulus of a complex number. In this paper, when we state that $\log|c|\approx 0,c\in\mathbb{C},$ we mean that $|c|$ is of order 1. In tropical geometry the map $\log|\cdot|$ is referred to as an Archimedean valuation. For a general introduction to tropical geometry we refer to [Sha11, MS15].

The Newton polytope of $f$ is the line segment $[0,d]\subset\mathbb{R}$. The convex hull of the points $\{(i,v_{i})\}_{0\leq i\leq d}\subset\mathbb{R}^{2}$ is called the lifted Newton polytope. We will consider the upper hull of the lifted Newton polytope. For a specific example, this is shown as a solid black line in Figure 3. The vertices of this upper hull are the points $(\beta_{\ell},v_{\beta_{\ell}})$, $\ell=0,\ldots,s$, with

$$0=\beta_{0}<\beta_{1}<\ldots<\beta_{s}=d\,.$$

We call the set $\Delta=\{(\beta_{0},\beta_{1}),(\beta_{1},\beta_{2}),\ldots,(\beta_{s-1},\beta_{s})\}$ the subdivision induced by the coefficients $c_{i}$, or the induced subdivision for short. We say that a point $\tau\in\mathbb{R}\cup\{-\infty\}$ is a root of $\texttt{t}{f}$ if the maximum in (8) is attained at least twice. A root of $\texttt{t}{f}$ is called a tropical root of $f$. The multiplicity of a root $\tau$ of $\texttt{t}{f}$ is the number $\beta_{\ell}-\beta_{\ell-1}$ where $\beta_{\ell}$ and $\beta_{\ell-1}$ are the largest, respectively the smallest value of $i$ for which $v_{i}+i\tau$ is maximal. Counted with multiplicity, $\texttt{t}{f}$ has $d$ roots $\tau_{1},\ldots,\tau_{d}$ and we can give a closed formula for them. For $\beta_{\ell-1}<i\leq\beta_{\ell}$ we have

$$\tau_{i}=\frac{1}{m_{\ell}}(v_{\beta_{\ell-1}}-v_{\beta_{\ell}})=\log\;\left|\frac{c_{\beta_{\ell-1}}}{c_{\beta_{\ell}}}\right|^{\frac{1}{m_{\ell}}}$$

where $m_{\ell}=\beta_{\ell}-\beta_{\ell-1}$ is the multiplicity. In particular, the definition implies

$$\tau_{1}=\tau_{\beta_{0}}=\cdots=\tau_{\beta_{1}}<\tau_{\beta_{1}+1}=\cdots=\tau_{\beta_{2}}<\cdots<\tau_{\beta_{s-1}+1}=\cdots=\tau_{\beta_{s}}=\tau_{d}.$$

Tropical roots of polynomials are used for scaling (matrix) polynomial eigenvalue problems, see for instance [BNS13, GS09, NST15].

Furthermore, for $i=0,\ldots,d$ we define the constants

$$r_{i}=\begin{cases}1,&i=\beta_{\ell}\text{ for some }\ell\\ \exp(v_{\beta_{\ell}}+(\beta_{\ell}-i)\tau_{i}-v_{i}),&\beta_{\ell-1}<i<\beta_{\ell}\text{ for some }\ell\text{ and }v_{i}\in\mathbb{R}\\ \exp(v_{\beta_{\ell}}+(\beta_{\ell}-i)\tau_{i}),&\beta_{\ell-1}<i<\beta_{\ell}\text{ for some }\ell\text{ and }v_{i}=-\infty\end{cases}\,.$$

Geometrically, if $c_{i}\neq 0$ then $\log r_{i}$ is the distance of $v_{i}=\log|c_{i}|$ to the upper convex hull of the lifted Newton polytope. Figure 3 illustrates these concepts. Note that $\tau_{i}$, $\beta_{\ell-1}<i\leq\beta_{\ell}$, is the negative slope of the line connecting $(\beta_{\ell-1},v_{\beta_{\ell-1}})$ and $(\beta_{\ell},v_{\beta_{\ell}})$. We note that the complexity of computing the tropical roots of $f$ is linear in the degree $d$, see e.g. [Sha11, Proposition 2.2.1].

In the following we show that an element-wise mixed backward error of order machine precision also implies a tropical backward error of order machine precision and that the converse holds, under suitable assumptions, as well. As in Section 2 for the quadratic case, we assume in the following that $|x_{j}|$, $|\hat{x}_{j}|$ and $|\tilde{x}_{j}|$ have the same order of magnitude for $j=1,\ldots,d$.