Incorrect Left and Right Eigenvectors in Mathematica

Question

I'm trying to find the left and right eigenvectors of a pretty straightforward matrix, but Mathematica doesn't seem to be able to do it for even a 200 dimensional matrix.

---

Background:

Given an operator $L$, the right and left eigenvectors corresponding to the same eigenvalue $\lambda$ satisfy $$L \phi_R = \lambda \phi_R$$ $$\phi_L^T L = \lambda \phi_L^T$$ respectively. We can take the transpose of the second equation to find that the left eigenvector satisfies $$L^T \phi_L = \lambda \phi_L$$ Now, it can be shown that the left and right eigenvectors for the same eigenvalue cannot be orthogonal $$\phi_L^T \phi_R \neq 0$$ I have discretised a system of coupled first order equations and am numerically trying to find the eigenvalues and eigenvectors but for some reason Mathematica returns left and right eigenvectors that are orthogonal. I tried increasing the precision beyond Machine Precision but that didn't seem to help either.

---

Code:

Initialise some variables

T = 10;
NN = 100;
nn = NN - 1;
λ = 1;
g = 1;
dt = 2 T/(NN - 1);
t[i_] := -T + (i - 1) dt

Create the matrix with dimensions 2 nn. The matrix has $(\lambda t,-\lambda t)$ alternating on the diagonals, $(g,0)$ alternating on both the upper and lower off-diagonals, $-\frac{1}{2\,dt}$ on the upper second off diagonal, and $\frac{1}{2\,dt}$ on the lower second off diagonal.

Mat= 
DiagonalMatrix[Flatten[Table[{λ t[i], -λ t[i]}, {i, 1, nn}]]] + 
DiagonalMatrix[Flatten[{Table[{g, 0}, {i, 1, nn - 1}], g}], 1] + 
Transpose[DiagonalMatrix[Flatten[{Table[{g, 0}, {i, 1, nn - 1}], g}], 1]] + 
DiagonalMatrix[Flatten[Table[{-1/(2 dt), -1/(2 dt)}, {i, 1, nn - 1}]], 2] - 
Transpose[DiagonalMatrix[Flatten[Table[{-1/(2 dt), -1/(2 dt)}, {i, 1, nn - 1}]], 2]];

For periodic boundary conditions, however, I also need to add entries in top right and the bottom left

Mat[[1, 2 nn - 1]] = 1/(2 dt);
Mat[[2, 2 nn]] = 1/(2 dt);
Mat[[2 nn - 1, 1]] = -1/(2 dt);
Mat[[2 nn, 2]] = -1/(2 dt);

Define the transpose (for left eigenvectors):

MatT = Transpose[Mat];

Find eigenvalues and eigenvectors:

System = Eigensystem[N[Mat]];

SystemT = Eigensystem[N[MatT]];

Find the corresponding left and right eigenvectors corresponding to the same eigenvalue. For example, the 198th eigenvalue of the matrix and the 196th eigenvalue of the transpose are the same, so the corresponding eigenvectors should NOT be orthogonal, but I find that

System[[2, 198]].Conjugate[SystemT[[2, 196]]]

returns 0.

---

I'm not really sure why Mathematica is doing this. Is it just finding the eigensystem incorrectly? Do I need to increase the step size (dt) significantly more? Just not sure why it's giving such obviously incorrect eigenvectors.

Answer 1

After some investigation I can mostly explain the (correct) behavior.

First off, we set this up using machine doubles. After computing the eigensystems we reorder by eigenvalues. While there is the possibility of messing up conjugate pairs (wherein right and left eigenvalues get ordered differently due to small numeric discrepancies), I checked and this does not happen here.

T = 10;
NN = 100;
nn = NN - 1;
\[Lambda] = 1;
g = 1;
dt = 2 T/(NN - 1);
t[i_] := -T + (i - 1) dt;
Mat = DiagonalMatrix[
    Flatten[Table[{\[Lambda] t[i], -\[Lambda] t[i]}, {i, 1, nn}]]] + 
   DiagonalMatrix[Flatten[{Table[{g, 0}, {i, 1, nn - 1}], g}], 1] + 
   Transpose[
    DiagonalMatrix[Flatten[{Table[{g, 0}, {i, 1, nn - 1}], g}], 1]] + 
   DiagonalMatrix[
    Flatten[Table[{-1/(2 dt), -1/(2 dt)}, {i, 1, nn - 1}]], 2] - 
   Transpose[
    DiagonalMatrix[
     Flatten[Table[{-1/(2 dt), -1/(2 dt)}, {i, 1, nn - 1}]], 2]];
Mat[[1, 2 nn - 1]] = 1/(2 dt);
Mat[[2, 2 nn]] = 1/(2 dt);
Mat[[2 nn - 1, 1]] = -1/(2 dt);
Mat[[2 nn, 2]] = -1/(2 dt);
Mat = N[Mat];
MatT = Transpose[Mat];
{vals, vecs} = Eigensystem[Mat];
order = Ordering[vals];
valsO = vals[[order]];
vecsO = vecs[[order]];
{valsT, vecsT} = Eigensystem[MatT];
orderT = Ordering[valsT];
valsTO = valsT[[orderT]];
vecsTO = vecsT[[orderT]];

One can readily check that the eigenvector sets correctly diagonalize the matrix. I repeated most computations after numericizing to 100 digits and everything worked there as well. So we have very good reason to believe Eigensystem is performing up to standards.

The claim is that corresponding eigenvectors should have nonzero Hermitian inner products, at least under the assumption that there are no multiple eigenvalues and all eigenvalues are nonzero (both of which hold for this example). That claim, suitably modified, is correct: we just need to work with conjugate pairs. To see this, let mu and nu be two nonconjugate eigenvalues, with vecmuR a right eigenvector for mu and vecnuL a left eigenvector for nu. I'll use . to denote the Hermitian inner product (as distinct from Mathematica Dot), * on the right to denote conjugates, and a space to denote multiplication by a scalar.

(vecnuL.Mat).vecmuR == nu vecnuL.vecmuR ==
vecnuL.(Mat.vecmuR) = vecnuL.(mu* vecmuR) == mu* vecnuL.vecmuR

Since mu* != nu by supposition, the inner product factor is zero. The upshot is that a given left eigenvector vecnuL gives a zero inner product with every right eigenvector not corresponding to the conjugate eigenvalue. As the right eigenvectors comprise a basis, it must give a nonzero dot product with the lone remaining eigenvector.

On checking one finds that this indeed holds. When I work with 100 digits and use Chop to remove things smaller than 10^(-90) I get results consistent with machine arithmetic but with much higher resolution. The strange finding is that several of the nonzero inner products are quite small. For example, we see that the 39th eigenvalue in the sorted list is real valued.

Position[vecsO, vec_ /; FreeQ[vec, Complex], {1}, 
 Heads -> False]

(* Out[295]= {{39}, {160}} *)

Since it is conjugate to itself, the corresponding left eigenvector is (Hermitian) orthogonal to every right eigenvector except the 39th. Below is at machine precision.

Max[Abs[vecsO[[39]].Conjugate[Transpose[Delete[vecsTO, 39]]]]]

(* Out[333]= 1.66588268076*10^-15 *)

At 100 digits the max becomes suitably smaller, indeed, zero to accuracy of 100 digits.

Max[Abs[vecsO[[39]].Conjugate[Transpose[Delete[vecsTO, 39]]]]]

(* Out[356]= 0.*10^-100 *)

So it must have nonzero inner product with the corresponding 39th right eigenvector. Both are real-valued (no surprise, sine the eigenvalue is real), so this is just the usual Dot.

vecsO[[39]].vecsTO[[39]]

(* Out[296]= \
8.64747245485202872314331535006214592917434236091583764774560978230566\
958588080450843494119*10^-10 *)

Indeed it is nonzero. Small, but nonzero.

I get a consistent result when I remove vecsTO[[39]], and compute the null vector for the remaining matrix of unit length.

Here is where I am surprised. Usually inner products of left and right paired eigenvectors indicates some measure of conditioning, with smaller meaning more ill-conditioned. Yet SingularValueList[Mat] and LUDecomposition[Mat][[3]] both seem to indicate a well-conditioned matrix. One thing to note, however, is that these measure conditioning for different purposes, with the eigenvector measure giving an estimate of proximity to a multiple eigenvalue. So maybe small inner products of corresponding eigenvector pairs just indicates that small perturbations of the matrix will move their eigenvalues into a set with multiplicity.

The above paragraph takes me pretty much to the border of familiarity with numeric linear algebra. So I cannot say if it explains suitably the set of small nonzero inner products. Maybe.

--- edit ---

One thing I later realized is that it is really the condition number of the matrix of eigenvectors that matters here. And it is not so well conditioned (it's on the order of 10^10). Given that fact, it is not surprising that there would be considerable orthogonality defects appearing. Which is what the dot products above indicate. Ergo, they (the inner products) are no longer such a surprise.

--- end edit ---