24
$\begingroup$

Bug introduced after 5.0, in or before 8.0 and persisting through 12.0.


I notice in the following example that wrong smallest 2 eigenvalues are resulted if calculating from a sparse matrix. But it gives correct result if we

  • calculate the smallest 3,4,... eigenvalues from the sparse matrix
  • calculate any smallest eigenvalues from the corresponding normal matrix

I found many cases with this behavior. Why and any remedy?

In the code below, the 104×104 matrix in the minimal example I found is imported from Pastebin because it exceeds the length limit of a post.

Import["https://pastebin.com/raw/PpDfY3EQ", "Package"];
mysparsemat = mymat;
mymat = Normal[mysparsemat];
m = 2;
Reverse@First[Eigensystem[mymat, -m]]
Reverse@First[Eigensystem[mysparsemat, -m]]
m = 4;
Reverse@First[Eigensystem[mymat, -m]]
Reverse@First[Eigensystem[mysparsemat, -m]]

The result of the above code is

{-0.712477 + 1.02863*10^-16 I, 0.712477 + 1.46577*10^-16 I}

{0.712477 - 1.44294*10^-11 I, 0.712656 - 2.12258*10^-11 I}

{-0.712477 + 1.02863*10^-16 I, 0.712477 + 1.46577*10^-16 I, 
-0.712656 - 1.05578*10^-16 I, 0.712656 + 6.49144*10^-16 I}

{-0.712477 - 5.10777*10^-10 I, 0.712477 - 5.44863*10^-12 I, 
 0.712656 + 3.10198*10^-11 I, -0.712656 + 3.64677*10^-10 I}

So the resultant wrong 2nd smallest eigenvalue is actually the correct 3rd or 4th smallest eigenvalue. (The eigenvalues should be real and doubly degenerate in absolute value as expected from the original problem's nature).

Update
As you might have noticed in the comments below, a first guess accusing the Arnoldi algorithm is irrelevant. (And Matlab gives the correct result.)


A bug report has been filed for this in the Wolfram community.

$\endgroup$
  • $\begingroup$ Please do not use the bugs tag until your observations have been confirmed to be a bug. $\endgroup$ – J. M.'s ennui Mar 30 '18 at 13:43
  • $\begingroup$ @J. M. Nice edit! Thank you so much! $\endgroup$ – xiaohuamao Mar 30 '18 at 13:55
  • 3
    $\begingroup$ I think the degeneracy is trapping the Arnoldi algorithm. Look at what happens for e.g. Eigenvalues[sm, -2, Method -> {"Arnoldi", "Criteria" -> "Magnitude", "StartingVector" -> SparseArray[1 -> 1, 104]}] and Eigenvalues[sm, -2, Method -> {"Arnoldi", "Criteria" -> "Magnitude", "StartingVector" -> SparseArray[52 -> 1, 104]}]. $\endgroup$ – J. M.'s ennui Mar 30 '18 at 14:06
  • 2
    $\begingroup$ As heretical as this will sound: if you have access to MATLAB, try using eigs() on your matrix. At the very least, it's a way to be sure that it's the algorithm itself that is going flat, and not Mathematica's specific implementation (even tho both are using ARPACK under the hood). $\endgroup$ – J. M.'s ennui Mar 31 '18 at 1:11
  • 7
    $\begingroup$ OK, I asked a friend with MATLAB to run this example. eigs() is perfectly capable of returning the two smallest eigenvalues (in magnitude). I don't know what's up with Mathematica now. $\endgroup$ – J. M.'s ennui Mar 31 '18 at 8:45
7
$\begingroup$

It appears that the proximity of the eigenvalues causes Eigensystem with the default parameters to be inaccurate. This can be fixed by increasing the basis size to 30

Import["https://pastebin.com/raw/PpDfY3EQ", "Package"]; 
mysparsemat = mymat; 
mymat = Normal[mysparsemat]; 
m = 2; 
Reverse[First[Eigensystem[mymat, -m]]]
Reverse[First[Eigensystem[mysparsemat, -m, Method -> {"Arnoldi", "BasisSize" -> 30}]]]

which then gives the desired result

{-0.712477 + 1.02679*10^-16 I, 0.712477 + 2.53294*10^-16 I}

{-0.712477 + 3.34996*10^-16 I, 0.712477 + 1.11772*10^-15 I}
$\endgroup$
0
$\begingroup$

I found a similar issue in a non-numeric caluclation:

sol = C[1] Cos[\[Lambda] x] + C[2] Sin[\[Lambda] x] + C[3] Cosh[\[Lambda] x] + C[4] Sinh[\[Lambda] x];
row1 = CoefficientArrays[sol /. {x -> 0}, {C[1], C[2], C[3], C[4]}][[2]];
row2 = CoefficientArrays[D[sol, x] /. {x -> 0}, {C[1], C[2], C[3], C[4]}][[2]];
row3 = CoefficientArrays[D[sol, x, x] /. {x -> l}, {C[1], C[2], C[3], C[4]}][[2]];
row4 = CoefficientArrays[D[sol, x, x, x] /. {x -> l}, {C[1], C[2], C[3], C[4]}][[2]];

sparsematrix = {row1, row2, row3, row4};
Det[sparsematrix]
normalmatrix = Normal@sparsematrix;
Det[normalmatrix]// FullSimplify

It gives 0 and 2 \[Lambda]^6 (1 + Cos[l \[Lambda]] Cosh[l \[Lambda]]). The first result is unexpected. It took me a while to figure that out that Det[] does different things and fails when sparse vectors are collected as normalform-matrix. It can be fixed by Det[SparseArray[sparsematrix]]. I consider this a bug. Det[] apparently understands that the matrix is 4 by 4 (leaving out one row in the above example causes Det[] to complain about the matrix not being square). So one assumes that it is handled properly, but the result is just wrong.

$\endgroup$
  • $\begingroup$ I doubt this is related to the original issue in this thread. It should be reported separately. $\endgroup$ – Daniel Lichtblau Sep 29 '20 at 4:48
  • $\begingroup$ I see. Where can I file a bug report? $\endgroup$ – Rainer Glüge Sep 30 '20 at 9:40
  • $\begingroup$ Actually I can just file it. I think for external reports it might be support@wolfram.com. Also poking around the web site I found this: wolfram.com/support/contact/email/?topic=feedback $\endgroup$ – Daniel Lichtblau Sep 30 '20 at 17:28
  • $\begingroup$ Report now filed. $\endgroup$ – Daniel Lichtblau Sep 30 '20 at 17:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.