Polynomial for very large number of roots

Question

Update: Please also see this solution here provided by MattL.

I have the roots for a very large order polynomial (>100), and from those alone wish to recreate the polynomial and run into numerical challenges when using the poly function for doing this available in Matlab and Python (numpy).

I attempted to convolve the coefficients for the roots and still see numerical issues (although I am not sure why as the dynamic range of the interim results when convolving subsets of the roots is far below that allotted by the double precision floating point numbers I am using).

As a simple example, and showing where I begin to see issues: I start with a 59th order polynomial with unity coefficients: $1+x+x^2+\ldots x^{59}$

The roots of this will be the 60th roots of unity given as $e^{jn2\pi/60}$ with n from 1 to 59, which I also get from Python using np.roots which are then given as complex double precision floating point representation (real and imaginary are both double precision). I then used the roots to poly function within python and would expect to get the unity coefficients back, but instead get large errors in the middle of the polynomial (yes, consistent with where we would have the biggest sum of products).

I then used my own approach of iteratively convolving coefficients while watching for dynamic range issues. For example in the plot below, I first show the resulting polynomial if only the first half of the roots were used, with the magnitude in dB to show that the dynamic range is far below the 300 dB+ allotted by double precision floating point. The lower portion of this plot shows the result of convolving this polynomial returned by the first 30 roots with the similar polynomial returned by the second 30 roots with no improvement in error.

I went on to iteratively convolve each of the complex roots $z_n$, by convolving the first two (as $x-z_1$ and $x-z_2$), and then convolving the result of that with the next root (as $z^2-(z_1+z_2)+z_1 z_2$ and $x-z_3$), etc through all of the roots, as given in Python below:

def cz(roots):
    z = list(roots)
    result = [1]
    while z:
        result = np.convolve(roots, [1, -z.pop()])
    return result

This resulted in even more error:

My question is two fold: what is the source of this numerical error and are their improved mathematical techniques that I could use for converting large numbers of roots (I would like to be able to extend this to 100s of roots) to its polynomial form? I used the roots of unity in this simple example, but more generally I would like to be able to do this with any arbitrary roots on the complex plane (and can limit it to roots that are on or inside the unit circle).

Answer 1

As I understand your issue, you are computing the coefficients of the product of two polynomials like in this example:

$$(ux^2+vx+w) \times (a_0 + a_1 x + a_2 x^2 + ...)=wa_0+(a_0v+a_1w)x+\cdots \tag{1}$$

by convolving their coefficients:

$$[u,v,w] \star [a_0, a_1, a_2, \cdots ]=[wa_0,(a_0v+a_1w),(a_0u+a_1v+a_2w),\cdots]\tag{2}$$

and you want to understand the origin of the instabilities you are finding.

I am going to propose you to write the convolution operation as the application of a certain banded matrix to a column matrix, in this way:

$$\begin{pmatrix} w&&&&&&&\\ v&w&&&&&&\\ u&v&w&&&&&\\ &u&v&w&&&&\\ &&\ddots&\ddots&\ddots &&&\\ &&&u&v&w&&\\ &&&&u&v&w&\\ &&&&&u&v&w \end{pmatrix}\begin{pmatrix}a_0\\a_1\\ a_2\\ \vdots \\ \vdots \\ a_k \\ \vdots\end{pmatrix}\tag{3}$$

In this presentation, the convolution of two sequences, which is a symmetric operation, has been made asymmetrical. It is a way to single out one of the two polynomials and understand its intrinsic rôle, whatever the other polynomial.

A nice way to consider (3) is by viewing the succession of lines as a "moving coefficient window" sliding in front of the second "table of values".

Please note that the column vector

• has been turned upside down as is classical for convolution.

• can be given "trailing zeros" in order to have the same number of elements as the number of columns of the matrix. Moreover, ending by a certain number of zeros (at least equal to the degree of the first polynomial) is usual for getting an exact convolution (these considerations are classical when one examines the sizes of the convolution of two sequences).

In general, if $d$ is the degree of the first polynomial, the banded matrix has $d+1$ bands.

Now a very bad case that will be a kind of counterexample displaying an intrinsic origin of the problem of instability you want to solve:

Let us consider polynomial

$$p(x)=x^5+3x^4+x^3+x^2+x+1\tag{4}$$ with degree $d=5$, a "perturbated" version of your polynomial. The associated banded matrix with $d+1=6$ bands with format $15 \times 15$ has the huge condition number $c=1.33 \times 10^7$ (recall : the condition number is the number by which the error is multiplied, here 10 million times). For example if we work with digits truncated to $10^{-8}$, the second rank decimals will be corrupted...).

If I consider the case of a $20 \times 20$ matrix, with the same polynomial, it is worse: the condition number is as high as $2 \times 10^9$ !

Moreover, we are still far from the objective of considering polynomials with degrees of some hundreds...

Had we taken polynomial $x^5+2x^4+x^3+x^2+x+1$ with a same $15 \times 15$ matrix, we would have had a different condition number: "only" $10^4$...

The conclusion is that the observed fluctuations depend heavily on the involved polynomial(s) (I use the plural ; a similar analysis could be done on the second polynomial). This does not mean that the rounding errors due to thousands of multiplications/additions do not play a rôle, but a secondary rôle.

Here is the way we can obtain the condition number with Matlab (please note that Matlab does almost all the work : the construction of the banded convolution matrix is handled by a specific function "convmtx"):

n=15;              % matrix size
p=[1,3,1,1,1,1,1]; % polynomial coeff.
M=convmtx(p,n);    % conv. matrix
M=M(:,1:n);        % truncation (otherwise M is rectangular)
M=M';              % optional in fact
cond(M)            % condition number

Remark 1: The roots of polynomial $p(x)$ given by (4) are

$$0.7729 \pm 0.8832i, \ \ -0.4685 \pm 1.1729i, \ \ -1.2426, \ \ -0.3663$$

Remark 2: A very interesting article here.

Remark 3: In order to understand how the condition number is used, see here.

Answer 2

(too long for a comment)

What I mean to say is that, if you have

$$ P = \prod\limits_{k = 1}^n {\left( {x - r_{\,k} } \right)} $$ then you may divide the roots into batches for which the coefficients formula is known and then convolve between them and with the remaining coefficients. That is

• if (some) of the roots are integral you have better to group them into $$ P_{int} = \prod\limits_{k = 1}^{n_{\,int} } {\left( {x - p_{\,k} } \right)} $$ and convolve the coefficients as $$ \begin{array}{*{20}c} 0 & 1 & 2 & \cdots & {n_{\,int} - 1} & {n_{\,int} } \\ \hline {} & { - p_{\,1} } & 1 & {} & {} & {} \\ {p_{\,2} p_{\,1} } & { - p_{\,2} } & 0 & {} & {} & {} \\ \hline {p_{\,2} p_{\,1} } & { - p_{\,1} - p_{\,2} } & 1 & {} & {} & {} \\ \end{array} $$ and so on, that finally leads to Vieta's formulas, which in case you can apply directly.
  The resulting coefficients will be exact (within the digits limit).

• if (some) of the roots are rational you can proceed as above, or you can extract the lcm of the denominators $$ P_{rat} = \prod\limits_{k = 1}^{n_{\,rat} } {\left( {x - {{p_{\,k} } \over {q_{\,k} }}} \right)} = \prod\limits_{k = 1}^{n_{\,rat} } {\left( {x - {{p'_{\,k} } \over Q}} \right)} = {1 \over {Q^{n_{\,rat} } }}\prod\limits_{k = 1}^{n_{\,rat} } {\left( {Qx - p'_{\,k} } \right)} $$

• if (some) of the roots are in aritmetic progression, then you know the coefficient for that batch through the Falling Factorial and the Stirling N. of 1st kind $$ \eqalign{ & P_{ap} = \prod\limits_{k = 1}^{n_{\,ap} } {\left( {x - kr} \right)} = \prod\limits_{k = 0}^{n_{\,ap - 1} } {\left( {x - r - kr} \right)} = r^{n_{\,ap} } \prod\limits_{k = 0}^{n_{\,ap} - 1} {\left( {{{x - r} \over r} - k} \right)} = \cr & = r^{n_{\,ap} } \left( {{{x - r} \over r}} \right)^{\,\underline {\,n_{\,ap} \,} } = r^{n_{\,ap} } \sum\limits_{0\, \le \,k\, \le \,n_{\,ap} } {\left( { - 1} \right)^{\,n_{\,ap} - k} \left[ \matrix{ n_{\,ap} \cr k \cr} \right]\left( {{{x - r} \over r}} \right)^{\,k} } = \cr & = r^{n_{\,ap} } \sum\limits_{0\, \le \,k\, \le \,n_{\,ap} } {\sum\limits_{0\, \le \,j\, \le \,k} {\left( { - 1} \right)^{\,n_{\,ap} - j} \left[ \matrix{ n_{\,ap} \cr k \cr} \right] \left( \matrix{ k \hfill \cr j \hfill \cr} \right)\left( {{x \over r}} \right)^{\,j} } } \cr} $$

• if (some) of the roots are in geometric progression, then also you know the coefficient for that batch via the q-binomial theorem $$ \eqalign{ & P_{gp} = \prod\limits_{k = 0}^{n_{\,gp} - 1} {\left( {x - r^{\,k} } \right)} = x^{n_{\,gp} } \prod\limits_{k = 0}^{n_{\,gp} - 1} {\left( {1 - {{r^{\,k} } \over x}} \right)} = x^{n_{\,gp} } \left( {1/x;\;r} \right)_{n_{\,gp} } = \cr & = x^{n_{\,gp} } \sum\limits_{0\, \le \,k\, \le \,n_{\,gp} } {\left( \matrix{ n_{\,gp} \cr k \cr} \right)_{\;r} \left( { - {1 \over x}} \right)^{\,k} r^{\binom{k}{2} } } \cr} $$

• if (some) of the roots lie on a circle with center at the origin and radius $r$, and are uniformly distributed with angle $0, 1/n 2 \pi , 2/n 2 \pi, \cdots , (n-1)/n 2 \pi$ then $$ \eqalign{ & x^{\,n} - r^{\,n} = 0\quad \Leftrightarrow \quad x = re^{\,i{k \over n}2\pi } = r\omega _n ^{\,k} \quad \left| {\;k \in \left\{ {0,1, \cdots ,n - 1} \right\}} \right.\quad \Leftrightarrow \cr & \Leftrightarrow \quad x^{\,n} - r^{\,n} = \prod\limits_{k = 0}^{n - 1} {\left( {x - re^{\,i{k \over n}2\pi } } \right)} \cr} $$ and if the center of the circle is at $c$ $$ \prod\limits_{k = 0}^{n - 1} {\left( {x - \left( {c + re^{\,i{k \over n}2\pi } } \right)} \right)} = \left( {x - c} \right)^{\,n} - r^{\,n} $$

Answer 3

For the leading example, one has to remark that the numpy.roots function is good at finding all the roots fast, but fails to give them to full accuracy. So inserting a Newton step greatly improves the roots and the reconstructed polynomial.

## in IPython, otherwise use print commands for output
p = [1.0]*60
z = np.roots(p)
max(abs(np.polyval(p,z)))
#>>>  1.069574264146017e-12
## Apply a Newton step,
## further steps have random effects with minimal improvements
z = z - np.polyval(p,z)/np.polyval(np.polyder(p),z)
max(abs(np.polyval(p,z)))
#>>> 4.6250852608698854e-14
## Reconstruct the polynomial
p1 = np.poly(z)
max(abs(p1-1))
#>>> 0.002289390374985656
plt.plot(p,'-+', ms=2); plt.plot(p1,'-o', ms=3); plt.show()

This is with one Newton step. A second one reduces the error by half, a third one increases the error again.

This is from the start better than any of your attempts. One might say that any polynomial that faithfully represents the roots with their errors needs to deviate from the original one in a systematic fashion. It is unlikely that errors, like rounding and cancellation, during the polynomial reconstruction will produce the exact opposite effect to the errors in the roots.