6
$\begingroup$

I would like to compute the roots of a polynomial with exponentially small coefficients.

$$ \sum_{n=0}^N a_n \frac{z^n}{\sqrt{n!}} \tag{$\ast$}$$

where $a_n$ are Normal random variables with mean $0$ and variance $1$. For simplicity you can just say $a_n \equiv 1$ always.

When the decay is large the coefficients are getting very small. $\sqrt{n!} \approx (n/e)^{n/2}$ this is causing underflow on my computer.


Q Calling np.roots on $\ast$ for $N=2000$ results errors because the coefficients are too small?.

What are some possible ways to handle numbers that are getting small this quickly, or can the power method simplify in some other way in this case?


What I Learned

Numpy's np.roots appproximates roots as the eigenvalues of the companion matrix of the polynomial.

$$ \left[ \begin{array}{cccl} 0 & 0 & \dots & 0 & -a_n \\ 1 & 0 & \dots &0 & -a_{n-1} \\ 0 & 1 & \dots & 0 & -a_{n-2} \\ 0 & \vdots & & \vdots & \vdots \\ 0 & 0 & \dots & 1 & -a_0 \end{array}\right]$$

In term, these eigenvalues are computed using something called LAPACK.

Without referring to a library, I could maybe implement the power method on my own. From initial random vector $b_0 \in S^N$, multiply by $A$ and rescale.

$$ b_{k+1} = \frac{A b_k}{|| A b_k ||}$$

The coefficients are getting small very quickly but the ratio of neighbor coefficients is not growing too fast $b_n = a_n/\sqrt{n!}$

$$ b_n/b_{n+1} = \sqrt{n+1} $$


Really Simple Example

What are the roots of $\epsilon x^2 + ax + b = 0$ ? The quadratic formula returns two numbers approximately

$$ x = \frac{ -a + \sqrt{a^2 - 4b \epsilon} }{2\epsilon } = \frac{ -a + a\sqrt{1 - 4(b/a) \epsilon} }{2\epsilon } \approx -\frac{a}{\epsilon} , -\frac{b}{a}$$

The ratio of the coefficients seems to matter a lot. In our case, they are of the same order of magnitude $\frac{a}{\epsilon} , \frac{b}{a} \approx \sqrt{n}$


Here is some code I wrote incorporating the rescaling suggestion... not many lines. Certainly it's enough to get the roots, but not sure if I trust it for finer statistics. Convergence is king here, I think.

N = 2000
a = 0.5*(np.random.normal(0,1,N) + 1j*np.random.normal(0,1,N))

fact = np.cumproduct(np.sqrt(np.arange(N)+1)/np.sqrt(N/np.exp(1)))

w = np.roots(a/fact[::-1])
plt.plot(w.real, w.imag, '.')

enter image description here


Final Version

The obvious (in hindsight) rescaling leads to much more accurate and convincing pictures, that I believe will hold up to further scrutiny.

N = 2000
a = 0.5*(np.random.normal(0,1,N) + 1j*np.random.normal(0,1,N))

A = 1j*np.zeros((N,N))
A[:,-1] = a

x = np.arange(N-1)
A[x+1, x] =np.sqrt(x+1)

w = np.linalg.eigvals(A)
plt.plot(w.real, w.imag, '.')

enter image description here

$\endgroup$
  • $\begingroup$ Did you mean to say "roots of a polynomial with exponentially small coefficients" in the first sentence? $\endgroup$ – Bill Barth Feb 13 '15 at 14:26
  • 4
    $\begingroup$ $\sqrt{N!}$ is not representable in double precision arithmetic for $N=2000$ (Sterling's formula indicates that it is roughly equal to $(2000/e)^{1000}\approx 736^{1000}>10^{2800}$). But even if it were, any kind of floating point arithmetic would suffer from catastrophic loss of accuracy because you are evaluating and adding numbers that are on vastly different scales. $\endgroup$ – Wolfgang Bangerth Feb 13 '15 at 15:09
  • $\begingroup$ Have you tried the rescaling trick I suggested in a comment to your other question? That could help. $\endgroup$ – Federico Poloni Feb 13 '15 at 21:40
  • $\begingroup$ @FedericoPoloni The scaling constants are super-exponentially decreasing, so why would that solve the problem? $\endgroup$ – Kirill Feb 13 '15 at 22:09
  • $\begingroup$ @Kirill I am not hoping to eliminate completely the superexponential growth problem, just to mitigate it enough for the numerical methods to converge. $\endgroup$ – Federico Poloni Feb 13 '15 at 22:24
7
$\begingroup$

Note that, if $D$ is invertible, the eigenvalues of $A$ and $DAD^{-1}$ are the same.

You can avoid floating-point underflow when forming the matrix by scaling the companion matrix by a diagonal matrix $D$ such that $D_{ii} = 1/\sqrt{(n-i)!}$. For the polynomial $x^4 + a x^3/\sqrt{1!} + b x^2/\sqrt{2!} + c x/\sqrt{3!} + d/\sqrt{4!}$, this gives a modified companion matrix of

$$ \left(\begin{array}{cccc} & & & -d\\ 1/\sqrt{4} & & & -c\\ & 1/\sqrt{3} & & -b\\ & & 1/\sqrt{2} & -a\\ \end{array}\right). $$

You can use a trick such as $p(x) \to x^n p(1/x)$ that inverts the roots of $p$ with respect to the circle in order to get $p$ into the form given above.

I do not know whether LAPACK reliably finds the eigenvalues of such matrices, but all of the entries here are quite reasonably scaled.

$\endgroup$
  • $\begingroup$ Of course, for a companion matrix built from coefficients with widely varying magnitudes, a preliminary balancing is necessary. An alternative is to instead build a matrix pencil (thus avoiding the need to build a monic polynomial) and use the QZ algorithm with balancing. $\endgroup$ – J. M. Nov 18 '15 at 11:46

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.