45
$\begingroup$

Bug introduced in 12.0 and fixed in 12.1


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$
20
  • 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$ Commented May 18, 2019 at 0:50
  • 14
    $\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
    Commented May 18, 2019 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
    Commented May 18, 2019 at 4:05
  • 8
    $\begingroup$ Filing a report (apparently that has not been done to date). $\endgroup$ Commented Jun 10, 2019 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
    Commented Jul 29, 2019 at 14:04

4 Answers 4

10
$\begingroup$

Fixed in 12.1


Mathematica graphics


ClearAll[x, n];
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) // Chop
$\endgroup$
1
  • 2
    $\begingroup$ Please edit and update the bug header too. $\endgroup$
    – Szabolcs
    Commented Mar 18, 2020 at 18:19
32
$\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$
1
  • 3
    $\begingroup$ This persists for Eigensystem also. $\endgroup$ Commented May 19, 2019 at 3:38
12
$\begingroup$

Edit: Eigensystem fixed in 12.1 in addition to Eigenvalues

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 (?!))
(We waited 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$
10
  • $\begingroup$ The generated eigenvectors aren't even really eigenvectors of the matrix, so the error is not just in the eigenvalues. $\endgroup$
    – Roman
    Commented May 19, 2019 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$ Commented May 19, 2019 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
    Commented May 19, 2019 at 19:35
  • 6
    $\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
    Commented May 19, 2019 at 20:44
  • 13
    $\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$ Commented May 22, 2019 at 2:29
4
$\begingroup$

I am not qualified to be at this site because the last time I used eigenvectors was well over half a century ago. The word "stiff matrix" came back to me, so I increased the precision of the author's code by rounding the two real numbers to 50 places. It took forever to compute, but Mathematica solved the problem accurately. That is, R = N[5/100, 50]; and mat0 = mat[N[-2/10 [Pi], 50]]; I ended up with zero to 47 places.

Occasionally--when solving differential equations numerically--I came across stiff systems, so I checked for this condition prior to working with them (I forget the method I used).

Again, sorry for my layman's interjection.

Nick Bagley

$\endgroup$
1
  • 4
    $\begingroup$ Charles, according to Wikipedia, "[...] the eigenvalue problem for all normal matrices is well-conditioned." So stiffness should not play a role here, unless a bad algorithm is being used. $\endgroup$
    – Roman
    Commented Jan 26, 2020 at 17:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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