mathnet / mathnet-numerics

Math.NET Numerics
http://numerics.mathdotnet.com
MIT License
3.51k stars 898 forks source link

Eigenvalue Calculation does not finish #1090

Open sfetzel opened 3 months ago

sfetzel commented 3 months ago

Consider the following code, which calculates the roots of the polynomial x^4 using the companion matrix eigenvalues:

var polynomial = new Polynomial(0, 0, 0, 0, 1);
var rootsReal = polynomial.EigenvalueMatrix().Evd(Symmetricity.Asymmetric).EigenValues.ToArray();
var rootsComplex = polynomial.EigenvalueMatrix().ToComplex().Evd(Symmetricity.Asymmetric).EigenValues.ToArray();

The second line calculates the roots using the eigenvalues of the real companion matrix and the last line uses the complex companion matrix.

Actual behavior The execution of the last line of code will not finish.

Expected behavior The execution should termiante and the contents of the variable rootsComplex should be equal to the contents of the variable rootsReal.

anshuman0se commented 1 month ago

Hi, I can work on this

anshuman0se commented 1 month ago

Hello, after further inspection, I believe the issue is to do with how complex matrices are treated by the EVD algorithm in this library. When matrices are converted to complex form, certain properties may no longer hold, so many issues can occur (such as infinite loops occuring when calculating eigenvalues).

In your code you might not handle matrix symmetry correctly. This could alter the results of the calculations performed by the algorithm. You can check if the matrix is symmetric to see which EVD algorithm to use. Here is my corrected version of your code:

var polynomial = new Polynomial(0, 0, 0, 0, 1);
var rootsReal = polynomial.EigenvalueMatrix().Evd(Symmetricity.Asymmetric).EigenValues.ToArray();
// Convert the companion matrix to a complex matrix.
var complexMatrix = polynomial.EigenvalueMatrix().ToComplex();

// Check if the complex companion matrix is symmetric.
if (complexMatrix.IsSymmetric())
{
// Use symmetric eigenvalue decomposition if the matrix is symmetric.
var rootsComplex = complexMatrix.Evd(Symmetricity.Symmetric).EigenValues.ToArray();
}
else
{
// Use asymmetric eigenvalue decomposition if the matrix is not symmetric.
var rootsComplex = complexMatrix.Evd(Symmetricity.Asymmetric).EigenValues.ToArray();
}
sfetzel commented 1 month ago

Thank you for your response.

Afaik the result of polynomial.EigenvalueMatrix is the companion matrix which is asymmetric. I also specified this correctly in my code. The matrix looks like this:

<-0; 0>  <-0; 0>  <-0; 0>  <-0; 0>
 <1; 0>   <0; 0>   <0; 0>   <0; 0>
 <0; 0>   <1; 0>   <0; 0>   <0; 0>
 <0; 0>   <0; 0>   <1; 0>   <0; 0>

The code you wrote does not terminate as well.

Please note, that this does not happen when the Intel MKL provider is used (I did not test OpenBLAS).

Afaik eigenvalue calculations can get stuck. Implementations must be aware of this and add a technique to get rid of this.

MH-ZShearer commented 1 month ago

Sorry, I misunderstood how the providers worked for Math.NET so the information I provided was not relevant. I am still interested in a resolution to this though.

We are currently experiencing issues trying to factor the roots of some polynomials. The MKL and other libraries can accomplish it correctly, but we can't justify the 150+ MB MKL package just for one method.