36
$\begingroup$

Bug introduced in 12.0.


The following code calculates the eigenvalues of a certain complex matrix, which come in pairs of opposite complex numbers. Therefore one can check whether the sum of all eigenvalues is equal to the trace of the matrix, which is zero.

This is indeed the case in Version 10.1 & 11.3 as far as I tested. However, Version 12.0 (Windows, Mac, Linux) gives something seriously wrong.

NN = 374; R = 0.05;
t1 = -1 + Cos[x] - I Sin[x] + I R; t1p = -1 + Cos[x] + I Sin[x] + 
  I R;
mat[x_] = 
  DiagonalMatrix[Table[If[EvenQ[n], t1, -1], {n, 0, 2 NN - 1 - 1}], 
    1] + DiagonalMatrix[
    Table[If[EvenQ[n], t1p, -1], {n, 0, 2 NN - 1 - 1}], -1] + 
   DiagonalMatrix[Table[If[EvenQ[n], -1, 0], {n, 0, 2 NN - 1 - 3}], 
    3] + DiagonalMatrix[
    Table[If[EvenQ[n], -1, 0], {n, 0, 2 NN - 1 - 3}], -3];
mat0 = mat[-0.2 \[Pi]];
Tr@mat0  (* 0. *)
Total@Eigenvalues@mat0  (* 0.394003 - 0.566499 I *)

I would rather switch back to 11.3 for a while. This looks really dangerous...


Original post of a more complex matrix with the same issue:

The code plots the real part of adding each pair. So the correct plot should be just zeros everywhere. This is the case in Version 10.1 & 11.3 as far as I tested (scattered numbers around $10^{-14}$ or so). However, Version 12.0 (Windows, Mac, Linux) gives something different as shown below.

NN = 200; R = 0.05;
xlist = Table[x, {x, -0.2 \[Pi], 0.2 \[Pi], 0.01}];
modl[n_] := 2*^-3 (Quotient[n, 2] - NN/2);
t1 = -1 + Cos[x] - I Sin[x] + I R; t1p = -1 + Cos[x] + I Sin[x] + I R;
t2a[n_] := -1 - modl[n]; t2b[n_] := -1 + modl[n];
mat[x_] = 
  DiagonalMatrix[
    Table[If[EvenQ[n], t1, t2a[n]], {n, 0, 2 NN - 1 - 1}], 1] + 
   DiagonalMatrix[
    Table[If[EvenQ[n], t1p, t2a[n]], {n, 0, 2 NN - 1 - 1}], -1] + 
   DiagonalMatrix[
    Table[If[EvenQ[n], t2b[n], 0], {n, 0, 2 NN - 1 - 3}], 3] + 
   DiagonalMatrix[
    Table[If[EvenQ[n], t2b[n], 0], {n, 0, 2 NN - 1 - 3}], -3];
list0 = Sort@Re@Eigenvalues[mat[xlist[[3]]]];
list0p = Table[list0[[i]] + list0[[2 NN - i + 1]], {i, NN}];
ListPlot[Tooltip@list0p, PlotRange -> All]

enter image description here

$\endgroup$
  • 1
    $\begingroup$ bugs: "This tag is reserved for questions where the problem has been vetted by this community and the observed behavior is confirmed to be a bug. Please do not use this tag for new questions." $\endgroup$ – AccidentalFourierTransform May 18 at 0:50
  • 11
    $\begingroup$ A simpler test is comparing Tr[matrix] with Total[Eigenvalues[matrix]] which shows a large discrepancy. I think this is worth reporting to support. It happens on OSX as well. $\endgroup$ – Carl Woll May 18 at 1:14
  • 1
    $\begingroup$ @Carl Woll. Interesting: On my machine (Windows 10; version 12.0) Tr[matrix] gives 0. and Total[Eigenvalues[matrix]] gives -2.23821*10^-13 - 2.66454*10^-15 I. $\endgroup$ – ulvi May 18 at 4:05
  • 4
    $\begingroup$ Filing a report (apparently that has not been done to date). $\endgroup$ – Daniel Lichtblau Jun 10 at 23:14
  • 1
    $\begingroup$ @xiaohuamao To my understanding, the problem is with the MKL library. It has been fixed in the development version (which uses a newer version of MKL). $\endgroup$ – ilian Jul 29 at 14:04
25
$\begingroup$

Not a solution but too big for a comment. There seems to be a catastrophic failure in Eigenvalues happening that is not due to the matrix being crazy. As a diagnostic, let's calculate the smallest (by absolute value) eigenvalue of the upper-left $n\times n$ part of the matrix

M = mat[xlist[[3]]];

For odd $n$ the answer is zero, so let's only do this for even $n$. We do this in two ways

  1. Calculate all eigenvalues and pick the one with the smallest absolute value:
    e1[n_?EvenQ] := M[[;; n, ;; n]] // Eigenvalues // Abs // Min
  1. Calculate only the smallest eigenvalue (by absolute value) with the Arnoldi algorithm:
    e2[n_?EvenQ] := Eigenvalues[M[[;; n, ;; n]], 1, 
      Method -> {"Arnoldi", "Criteria" -> "Magnitude", "Shift" -> 0}] // First // Abs

Method (2) is very reliable, whereas method (1) breaks down for $n=358$ and above:

enter image description here

Considering that the Arnoldi algorithm has no problems with this matrix, there seems to be something really strange going on in method 1.

$Version
(* 12.0.0 for Mac OS X x86 (64-bit) (April 7, 2019) *)
$\endgroup$
  • $\begingroup$ This persists for Eigensystem also. $\endgroup$ – CA Trevillian May 19 at 3:38
7
$\begingroup$

I attempted a workaround, to see if Eigensystem had any issues also. It does. This is very unfortunate.

(Will we have to wait for 12.1 for the fix (?!))

My code here:

e3[n_?EvenQ] := Eigensystem[M[[;; n, ;; n]]][[1]] // Abs // Min

Produces the following, which matches with @Roman shows:

Issue persisting for Eigensystem

(Apologies the colors/styles don't match with the plot from @Roman !!)

$Version
(* 12.0.0 for Microsoft Windows (64-bit) (April 6, 2019) *)
$\endgroup$
  • $\begingroup$ The generated eigenvectors aren't even really eigenvectors of the matrix, so the error is not just in the eigenvalues. $\endgroup$ – Roman May 19 at 10:32
  • $\begingroup$ If I understand what you mean, this would be some fundamental (internal) issue with the eigensolver being used? Maybe some error in the linear decomposition being used? $\endgroup$ – CA Trevillian May 19 at 17:48
  • $\begingroup$ Yes it looks like a serious foul in the default "Direct" diagonalizer. This is very surprising for something that's been around for such a long time (coded in LAPACK). $\endgroup$ – Roman May 19 at 19:35
  • 5
    $\begingroup$ No need to fear, there are easy sanity checks, as @CarlWoll has already posted one (trace rule). A sanity check for the eigenvectors is to compute {eval, evec} = Eigensystem[M] and then go evec.Transpose[M].Inverse[evec] - DiagonalMatrix[eval] // Norm which should give zero (or something close to it) if the eigenvectors evec truly diagonalize the matrix M. $\endgroup$ – Roman May 19 at 20:44
  • 10
    $\begingroup$ This sort of breakage of previously working code now happens so frequently on major updates that I've decided to keep one or more old executables of Mathematica on hand. Note: test suites should be catching all of these! $\endgroup$ – Ralph Dratman May 22 at 2:29

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.