19
$\begingroup$

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).

np.poly

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.

poly through convolution

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 version of conversion

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).

$\endgroup$
11
  • 2
    $\begingroup$ The most common numerical issue is the Runge-phenomenon , the polynomial oscillates extremely heavily. But this rather occurs if we have many equidistant real roots. No idea what went wrong in the case of the polynomials you chose. $\endgroup$
    – Peter
    Commented Dec 4, 2021 at 15:29
  • 3
    $\begingroup$ Do you compute $\mathrm{e}^{\mathrm{j} n 2\pi/60}$ by raising $\zeta_{60} = \mathrm{e}^{\mathrm{j} 2\pi/60}$ to various integer powers (so all of your roots have systematic copies of any (rounding/truncation) error in $\zeta_{60}$) or do you calculate $\mathrm{j} n 2\pi/60$ for each root and then exponentiate that? It might be worthwhile to look at how many steps of random walk of step-size "error in $\zeta$" the coefficient of $x^{30}$ takes in the giant sum of products of various powers of $\zeta$. (Related: never do sensitivity analysis on a discrete Fourier transform...) $\endgroup$ Commented Dec 4, 2021 at 15:52
  • $\begingroup$ @Peter Interesting, I will need to learn more about the Runge-phenomenon, as "oscillating" feels like what I am seeing. My roots are equidistant on the unit circle so may be a similar effect. Thanks for sharing that. $\endgroup$ Commented Dec 4, 2021 at 15:59
  • 1
    $\begingroup$ @JeanMarie No I haven't and what you are saying sounds very promising but is beyond my comprehension. Could you provide that idea as an answer here with a little more detail that I can then try? (I'm an engineer, not a mathematician so you are starting to venture into a foreign language for me- but even if you point me in the direction of what I do need to learn in order to be able to solve something like this I will be very appreciative) $\endgroup$ Commented Dec 4, 2021 at 16:52
  • 2
    $\begingroup$ All right. I will do that. For the companion matrix definition see here. $\endgroup$
    – Jean Marie
    Commented Dec 4, 2021 at 16:59

3 Answers 3

10
$\begingroup$

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.

$\endgroup$
8
  • $\begingroup$ Thank you so much! I will try this and comment on how it worked for me (I likely won’t get back to this until tomorrow) $\endgroup$ Commented Dec 4, 2021 at 23:26
  • $\begingroup$ +1, interesting. I've always seen the condition number used to estimate the error in $x$ in the equation $Ax = b$ given an error in $b$. This is what is done in the MathWorks link in Remark 3 and in the Wiki article on condition numbers. That is, a very large condition number means a relatively small error in $b$ can become a relatively large error in $x$. However, I feel in this case you are estimating the error in $b$ given the error in $x$. I'm not sure condition numbers give a good estimate; shouldn't the error in $b$ be bounded by the largest singular value of $A$ times the error in $x$? $\endgroup$
    – user196574
    Commented Dec 5, 2021 at 12:32
  • $\begingroup$ @user196574 Your remark about the "knowing b, find x" is interesting because it says implicitly that the "true usage" of condition number here would be to quantitfy the error in a "deconvolution" operation... I must think about that. My idea, you have understood it, is to isolate the contribution of a polynomial. The condition number comes indeed from considerations about linear system's solving, but has been recognized as general indicator of "good behavior" for a matrix. But for sure I don't pretend there are not better indicators... $\endgroup$
    – Jean Marie
    Commented Dec 5, 2021 at 12:42
  • 1
    $\begingroup$ The more I think on it, the more I think my comment was confused. The fact that the condition number is the same for $A$ and $A^{-1}$ means that there is strong symmetry between "given error in $b$, find error in $x$" and "given error in $x$, find error in $b$". Also, in my comment above I said that "error in $b$" is bounded by the largest singular value of $A$ times the error in $x$; that's just absolute error: the relative error will be largest if $A$ shrinks the magnitude of $x$ and stretches $\delta x$, so the relative error is determined by the ratio of largest and smallest sing values. $\endgroup$
    – user196574
    Commented Dec 5, 2021 at 12:46
  • 1
    $\begingroup$ So to summarize, although I was first hesitant on the condition number being used in this way, I now am pretty happy with it as providing a bound for the relative error. I expect that when used in the case $Ax = b$ to estimate the relative error in $b$ given the relative error in $x$, the actual relative error will be close to the upper bound provided by the condition number when $Ax$ has smaller magnitude than $x$. Here I am assuming that the large condition numbers are stemming from small singular values, so that the worst relative error occurs when $x$ is shrunk by those small sing values. $\endgroup$
    – user196574
    Commented Dec 5, 2021 at 12:59
4
$\begingroup$

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()

enter image description here

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.

$\endgroup$
1
  • $\begingroup$ Oh that’s a really good find as well thank you! This makes a lot of sense. @MattL $\endgroup$ Commented Dec 5, 2021 at 12:21
4
$\begingroup$

(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} $$

$\endgroup$
4
  • $\begingroup$ Thank you! This is interesting. All of my roots in the example given are complex as the roots of unity on the unit circle, so they all have magnitude 1 and angles given as $2\pi/60$ (in my example). That said I think they all are irrational. They are all however in complex conjugate pairs, so I could possibly group them into real 2nd order sections, (z^2-2z cos(w)+1), but I was also hoping for a more general solution (if it exists) as I will typically have irrational roots anywhere inside or on the unit circle. $\endgroup$ Commented Dec 4, 2021 at 19:39
  • 1
    $\begingroup$ I have added for the case of roots being the roots of unity ( or $r^n$). The most general formulas are the Vieta's formulas $\endgroup$
    – G Cab
    Commented Dec 5, 2021 at 15:25
  • $\begingroup$ @G Cab Hi ! Interesting. I have been interested by the reference (MattL) given in the front "Update" of the OP. This reference mentions the Leja orderin, I didn't know ; I haven't found the original paper but a paper that mentions it interesting on its own $\endgroup$
    – Jean Marie
    Commented Dec 10, 2021 at 18:36
  • $\begingroup$ @JeanMarie: thanks, that reference work is really interesting $\endgroup$
    – G Cab
    Commented Dec 10, 2021 at 19:29

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .