roots of polynomials of high degree: LinAlgError: Eigenvalues did not converge

Question

I wrote a simple script to generate random polynoimals $\displaystyle f(z)= \sum_{k=0}^N a_k \frac{z^k}{\sqrt{k!}} $ of high degree and find their roots. For more discussion on random polyomials see here.

N = 500

a = np.random.normal(0,1,N)

n = np.arange(N-1)+1
n = np.insert(n,0,1)
n = np.sqrt(n)
n = np.cumproduct(n)

z = np.roots(a/n[::-1])/np.sqrt(N)

plt.xlim([-1,1])
plt.ylim([-1,1])
plt.plot(z.real, z.imag, '.')

For smaller values of (here N=150) we do get a the picture of the unit disk filled with zeros:

For larger value of N there is an error message. I am sure NumPy can compute eigenvalues of larger matrices than this. Right?

LinAlgError: Eigenvalues did not converge

Answer 1

Using the inversion operation it is possible to circumvent some of the exponential growth issues of the coefficients:

Writing the coefficients backwards n[::1] instead of n[::-1] we get the reciprocals of the roots.

N = 500

a = np.random.normal(0,1,N)

n = (1.0 + np.arange(N-1))/N
n = np.insert(n,0,1)

n = np.sqrt(n)
n = 1.0/n
n = np.cumproduct(n)

z = np.roots(a*n[::1])
z = 1/z
L = 1.5
plt.xlim([-L,L])
plt.ylim([-L,L])
plt.plot(z.real, z.imag, '.')

t = np.arange(0,1,0.01)
z = 1*np.exp(2j*np.pi*t)
plt.plot(z.real, z.imag, '-')

Then using the function $\displaystyle z \mapsto \frac{1}{z}$ we can flip the roots back. Once the degree N is too large, this trick no longer works.