roots of polynomials with small coefficients

Question

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

---

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

Answer 1

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.