:

## Use Eigenvalues, not fsolve@Determinant

Maple 18

There has been a spate of Questions posted in the past week about computing eigenvalues. Invariably, the Questioners have computed some eigenvalues by applying fsolve to a characteristic polynomial obtained from a floating-point matrix via LinearAlgebra:-Determinant. They are then surprised when various tests show that these eigenvalues are not correct. In the following worksheet, I show that the eigenvalues computed by the fsolve@Determinant method (when applied to a floating-point matrix) are 100% garbage for dense matrices larger than about Digits x Digits. The reason for this is that computing the determinant introduces too much round-off error into the coefficients of the characteristic polynomial. The best way to compute the eigenvalues is to use LinearAlgebra:-Eigenvalues or LinearAlgebra:-Eigenvectors. Furthermore, very accurate results can be obtained without increasing Digits.

 The correct and incorrect ways to compute floating-point eigenvalues Carl Love 2016-Jan-18 restart: Digits:= 15: macro(LA= LinearAlgebra): n:= 2^5:  #Try also 2^3 and 2^4. A:= LA:-RandomMatrix(n): A is an exact matrix of integers; Af is its floating-point counterpart. Af:= Matrix(A, datatype= float): P:= LA:-CharacteristicPolynomial(A, x): P is the exact characteristic polynomial with integer coefficients; Pf is the floating-point characteristic polynomial computed by the determinant method. Pf:= LA:-Determinant(Af - LA:-DiagonalMatrix([x\$n])): RP:= [fsolve(P, complex)]: RP is the list of floating-point eigenvalues computed from the exact polynomial; RPf is the list of eigenvalues computed from Pf. RPf:= [fsolve(Pf, complex)]: RootPlot:= (R::list(complexcons))->      plot(           [Re,Im]~(R), style= point, symbol= cross, symbolsize= 24,           axes= box, color= red, labels= [Re,Im], args[2..]      ) : RootPlot(RP); RootPlot(RPf); We see that the eigenvalues computed from the determinant are completely garbage. The characteristic polynomial might as well have been x^n - a^n for some positive real number a > 1.   Ef is the eigenvalues computed from the floating-point matrix Af using the Eigenvalues command. Ef:= convert(LA:-Eigenvalues(Af), list): RootPlot(Ef, color= blue); We see that this eigenvalue plot is visually indistinguishable from that produced from the exact polynomial. This is even more obvious if I plot them together: plots:-display([RootPlot(Ef, color= blue), RootPlot(RP)]); Indeed, we can compare the two lists of  eigenvalues and show that the maximum difference is exceedingly small.   The following procedure is a novel way of sorting a list of complex numbers so that it can be compared to another list of almost-equal complex numbers. RootSort:= (R::list(complexcons))-> sort(R, key= abs*map2(`@`, signum+2, Re+Im)): max(abs~(RootSort(RP) -~ RootSort(Ef)));   