Accuracy of poly for random unitary matrices

I have been trying to do some computations with the characteristic polynomials of random unitary matrices, and I have come across some unusual behavior of the poly command.

Let's say that we generate a random complex unitary matrix and compute the coefficients of it's characteristic polynomial:

    n=400;
    [Q,R]=qr((randn(n)+1j*randn(n))/sqrt(2));
    U=Q*diag(diag(R)./abs(diag(R)));
    p=poly(U);
    plot(abs(p),'o');

Now, I believe that the polynomial coefficients should all be relatively close to one. However, the coefficients have magnitudes around 10^63 in some cases. This seems to be a problem for large n.

Now my question is if Matlab should compute the coefficients more accurately than it does? I would think that this problem is a relatively well conditioned one: recover the characteristic polynomial of a unitary matrix.

[José-Luis]

The joys of numerical precision

The larger n, the larger the chances that your data spans more orders of magnitude. That is bound to give matrix operations hiccups.

[John D'Errico]

You should not be surprised at this. Floating point arithmetic is NOT mathematics, although there are some occasional overlaps.

Answer 1

[Piers Lawrence]

It appears that the Matlab algorithm to compute the coefficients suffers from extreme cancellation in their formation. The magnitudes of the coefficients do not vary so much, and depending on the ordering in which they are formed, the intermediate values do vary a lot.

All of this is explained in On the Evaluation of Polynomial Coefficients which appeared in:

"D. Calvetti and L. Reichel On the evaluation of polynomial coefficients Numerical Algorithms, 33 (2003), pp. 153-161." 10.1023/A:1025555803588

The solution is to order the eigenvalues using leja ordering. The above code can be stabilized via:

    n=400;
    [Q,R]=qr((randn(n)+1j*randn(n))/sqrt(2));
    U=Q*diag(diag(R)./abs(diag(R)));
    lambda=eig(U);
    lambda_leja=leja(lambda);
    p=poly(lambda);
    p_leja=poly(lambda_leja);
    subplot(2,1,1);
    plot(abs(p),'o');
    title('Magnitude of polynomial coefficients (Unordered)')
    subplot(2,1,2);
    plot(abs(p_leja),'o')
    title('Magnitude of polynomial coefficients (Leja ordering)')