Which algorithm does MATLAB eig() use to diagonalize a complex symmetric matrix?

I used MATLAB eig() to find eigenvectors and eigenvalues of a complex symmetric matrix. I searched through MATLAB online documentation to find a link to the algorithm they use, but failed. Can someone link me to the algorithm used by MATLAB? My curiosity is piqued also because of the fact that the algorithm used by eig() didn't seem to be something simple enough. I am saying this because we have a rudimentary conjugate gradient complex symmetric eigensolver in FORTRAN, and we get poor quality of complex orthogonality* between eigenvectors, unlike MATLAB.

*note that for a complex symmetric matrix, eigenvectors corresponding to distinct eigenvalues have a zero transpose inner product, not a zero conjugate-transpose inner product. That is, $v_1^T v2=0$ , but $v†1v2≠0$.

Answer 1

[Stephan]

Hi,

Maybe you are looking for this:

https://de.mathworks.com/help/releases/R2018a/matlab/ref/eig.html?s_tid=doc_ta#btgapg5-1-algorithm

There are 2 algorithms used, depending on input and/or user choice in options.

Best regards

Stephan

[Pranav Gupta]

Sorry for the delay. I was unable to get my license number for a while, because of which I couldn't open the link you sent. I saw the link, but it doesn't talk about the specifics of the algorithm at all. Should I assume that the default algorithm is a QR algorithm? Again, I apologize for the late response.

[Bruno Luong]

According to Christine (who is TMW staff); it is LAPACK so more like Hessenberg reduction

LAPACK decision tree

Answer 2

[Christine Tobler]

EIG uses LAPACK functions for all cases. If there is a special case treatment for complex symmetric, I'm not aware of this.

Unless there are multiple eigenvalues, wouldn't a general nonsymmetric eigenvalue solver find eigenvectors that have a zero transpose inner product? (I haven't tried what EIG does for a complex symmetric matrix with multiple eigenvalues, because I'm not sure how to construct one).

[Bruno Luong]

You are right Christine, there is no warranty of transposed orthogonal of eigen vectors output for multiple eigenvalues case:

% Construct A complex symmetric with multiple eigen values
B=randn(5)+1i*rand(5);
[W,D]=eig(B.'*B);
A=W*diag([1 1 1 2 2]+1i*[2 2 2 3 3])*inv(W);
A = A + A.'; % make sure A is symmetric
disp(A)

   3.2186 + 5.0894i  -0.7773 - 0.7697i  -0.2351 - 0.1350i   0.2788 - 0.2281i   0.6007 + 0.6893i
  -0.7773 - 0.7697i   2.4414 + 4.5526i   0.1556 + 0.2721i  -0.0407 - 0.1935i  -0.2904 - 0.3498i
  -0.2351 - 0.1350i   0.1556 + 0.2721i   2.6024 + 4.1188i  -1.0149 - 0.6676i   0.3581 - 0.2851i
   0.2788 - 0.2281i  -0.0407 - 0.1935i  -1.0149 - 0.6676i   3.1806 + 6.1977i  -0.9264 - 0.1746i
   0.6007 + 0.6893i  -0.2904 - 0.3498i   0.3581 - 0.2851i  -0.9264 - 0.1746i   2.5570 + 4.0415i

[V,~]=eig(A);
V.'*V

ans =

   0.9406 + 0.0527i   0.0000 + 0.0000i  -0.0000 + 0.0000i   0.1090 - 0.0485i  -0.0000 + 0.0000i
   0.0000 + 0.0000i   0.9207 - 0.0703i   0.3896 + 0.0699i   0.0000 + 0.0000i   0.0186 - 0.0956i
  -0.0000 + 0.0000i   0.3896 + 0.0699i   0.7577 + 0.0086i  -0.0000 + 0.0000i  -0.1659 - 0.5657i
   0.1090 - 0.0485i   0.0000 + 0.0000i  -0.0000 + 0.0000i   0.7661 - 0.2008i   0.0000 - 0.0000i
  -0.0000 + 0.0000i   0.0186 - 0.0956i  -0.1659 - 0.5657i   0.0000 - 0.0000i   0.5822 - 0.3127i

[Christine Tobler]

Thanks, that's a good way to construct that matrix.