Eigenvectors of a general matrix

[lpawela]

As of #2445 eigenvalue calculation seems to work correctly. Yet, sometimes there is a problem with eigenvector computation. Consider the following examples.

let A = [[1, 2, 3], [2, 4, 0], [3, 0, 1]];
let cnt = 0.1;
let Ath = math.multiply( math.exp(math.multiply(math.complex(0,1), -cnt)), AA );
let Hth = math.divide(math.add(Ath, math.ctranspose(Ath)), 2);
math.eigs(Hth)

This code works correctly and produces correct eigenvalues. But if I change

let Hth = math.divide(math.add(Ath, math.transpose(Ath)), 2);
math.eigs(Hth)

I get the following error

Uncaught Error: Failed to find eigenvectors for the following eigenvalues: -2.32076242897952 + 0.2328529372998936i, 2.7657291798020682 - 0.27749853033260213i, 5.525058240845601 - 0.5543549068482564i
    at complexEigs.js:461:19
    at complexEigs.js:49:17
    at j (eigs.js:100:12)
    at Function.Array (eigs.js:51:14)
    at Object.eigs (typed-function.js:1085:85)
    at <anonymous>:1:6

The expected output would be

julia> A = [1 2 3; 2 4 0; 3 0 1];

julia> cnt = 0.1;

julia> Ath = exp(-1im*cnt)*A;

julia> Hth = (Ath + transpose(Ath))/2;

julia> eigen(Hth)
Eigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}
values:
3-element Vector{ComplexF64}:
 -2.320762428979524 + 0.23285293729989398im
 2.7657291798020984 - 0.27749853033260563im
  5.525058240845574 - 0.5543549068482575im
vectors:
3×3 Matrix{ComplexF64}:
  0.723537+0.0im           0.391404+4.34166e-17im  0.568592-7.08711e-17im
 -0.228519-3.1175e-17im   -0.641444-1.01409e-16im  0.732345+0.0im
 -0.651363-4.64016e-17im   0.659812+0.0im          0.374665+1.22787e-16im

The eigenvalues look correct, there is some problem with eigenvectors only. I tried some other examples and it seemed to work in that case. These were:

let A = [[3, -2], [math.complex(4, 2), -1]];
let Ah = math.ctranspose(A);
console.log(math.eigs(math.add(A, Ah)));
console.log(math.eigs(A));
console.log(math.eigs(Ah));

In this case there are no errors and the results are correct. I tried a handful of other Hermitian matrices and the error did not occur.

[josdejong]

Thanks for reporting. It would be nice to improve the math.eigs function further.

Anyone interested in looking into this? Help would be welcome

[zishiwu123]

I'm interested in looking into this. I have a question regarding the example given by @lpawela. What does the AA mean? I think it is a typo and you meant to have only one A.

let Ath = math.multiply( math.exp(math.multiply(math.complex(0,1), -cnt)), AA );

I started digging into src/function/matrix/eigs/complexEigs.js by adding some print statements to the function inverseIterate. Here is what I discovered from the first call to inverseIterate from findEigenvectors:

mathjs/src/function/matrix/eigs/complexEigs.js

Lines 445 to 446 in 99839a8

	while (solutions.length < multiplicities[i]) {
	const approxVec = inverseIterate(A, N, solutions, prec, type)

• orthog is an empty array
• vector b never gets an imaginary component
• norm(b) increases and decreases but never goes above 1.0 which is much smaller than the limit of largeNum=1000 so function returns null
  

  mathjs/src/function/matrix/eigs/complexEigs.js

  Lines 610 to 617 in 99839a8

  	let i = 0
  	while (true) {
  	b = randomOrthogonalVector(N, orthog, type)
  	b = usolve(A, b)
  	if (larger(norm(b), largeNum)) { break }
  	if (++i >= 5) { return null }
  	}

Here is the full print:

findEigenvectors: i = 0, multiplicities[i] = 1
inverseIterate
N: 3, orthog: , type: Complex
A
[
  [
    { re: 7.873676792142968, im: -0.790002779046487 },
    { re: -0.05449469718253663, im: 0.005467707572199387 },
    { re: -0.024544315085780977, im: 0.0024626458056940187 }
  ],
  [
    { re: -0.6633940925144037, im: 0.0665614287358509 },
    { re: 2.614862447423172, im: -0.2623613662107619 },
    { re: 0.015537913166105418, im: -0.0015589914224133727 }
  ],
  [
    { re: 70.87071457376122, im: -7.110789907219184 },
    { re: 515.887286883318, im: -51.76138176249111 },
    { re: 2.4437730390405754, im: -0.24519516652340045 }
  ]
]

i=0
b after randomOrthogonal
[
  { re: -0.5577373217901133, im: 0 },
  { re: -0.35390389747393836, im: 0 },
  { re: 0.750786994583116, im: 0 }
]
b after usolve
[
  [ { re: -0.07012143631526388, im: -0.007035611318852924 } ],
  [ { re: -0.13580167111455446, im: -0.013625616139934989 } ],
  [ { re: 0.30416250022046715, im: 0.030518044720311277 } ]
]
norm(b): 0.34211174277443424

i=1
b after randomOrthogonal
[
  { re: -0.48486519400867156, im: 0 },
  { re: -0.8741567859582652, im: 0 },
  { re: 0.027489219742454104, im: 0 }
]
b after usolve
[
  [ { re: -0.06322321085794975, im: -0.006343480129621954 } ],
  [ { re: -0.33103743476122005, im: -0.03321453246477565 } ],
  [ { re: 0.011136567183901904, im: 0.0011173838165543842 } ]
]
norm(b): 0.3388977341126987

i=2
b after randomOrthogonal
[
  { re: 0.9074931565986786, im: 0 },
  { re: 0.19886643261456535, im: 0 },
  { re: 0.3700112332155916, im: 0 }
]
b after usolve
[
  [ { re: 0.1150900999508549, im: 0.011547527438850812 } ],
  [ { re: 0.07440362437487062, im: 0.007465263253621657 } ],
  [ { re: 0.14990076095684626, im: 0.01504024369596464 } ]
]
norm(b): 0.2041253365902603

i=3
b after randomOrthogonal
[
  { re: 0.5595792179341783, im: 0 },
  { re: -0.7083943148664702, im: 0 },
  { re: -0.4301727484639601, im: 0 }
]
b after usolve
[
  [ { re: 0.06796888084805093, im: 0.006819635371904041 } ],
  [ { re: -0.2671750877955009, im: -0.02680692482336297 } ],
  [ { re: -0.1742736883344125, im: -0.017485693372155244 } ]
]
norm(b): 0.3277872524454639

i=4
b after randomOrthogonal
[
  { re: -0.14652413426547323, im: 0 },
  { re: -0.5243471598996081, im: 0 },
  { re: 0.8388031556825286, im: 0 }
]
b after usolve
[
  [ { re: -0.018752590363761725, im: -0.0018815350049009836 } ],
  [ { re: -0.20054639509272648, im: -0.02012175678954787 } ],
  [ { re: 0.3398200379947726, im: 0.03409573208027077 } ]
]
norm(b): 0.3970127437112941

[josdejong]

Thanks @zishiwu123 for looking into this 💪

[gwhitney]

It appears that some former author(s) may have neglected the fact that some of the auxiliary functions in complexEigs.js modify their leading argument by side effect, and/or lost track of the transforms that get the current working matrix back to the one originally supplied to the eigs function. Also, the only test case that might have found this had an incorrect reference value. The PR fixes these things and adds the example of @lpawela as a new test case. But I would say that the testing here is still a bit thin -- so if anyone looking at this issue has an additional test case or two that is going to exercise the non-symmetric eigenvector case, please file a PR adding it or at least let us know about the case. Thanks!