Accuracy of poly for random unitary matrices

Question

Piers Lawrence 2014-6-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/135173-accuracy-of-poly-for-random-unitary-matrices

回答： Piers Lawrence 2014-6-19

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.

2 个评论
显示无隐藏无

José-Luis 2014-6-18

编辑：José-Luis 2014-6-18

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 2014-6-18

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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Piers Lawrence 2014-6-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/135173-accuracy-of-poly-for-random-unitary-matrices#answer_141648

在 MATLAB Online 中打开

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