How to get eigenvalues of matrix (without eig, eigs) ?

I try to make a program, for identification of eigenvalues of matrix M (without eigs, eig, ... ).

I use QR algorithm, which should work for Hessenberg matrices.

(wiki: https://en.wikipedia.org/wiki/Eigenvalue_algorithm)

clc, clear;      % clear the window and workspace
M = magic(6);    % some sample matrix
ideal = eigs(M)  % get ideal solution by original matlab function
hM = hess(M);    % converting matrix to Hessenberg form
eigs(hM)  % checking, if eigenvalues of hM are the same.. as in M

iterations = 800;
A = hM;
% using QR algorithm for detection eigenvalues
for k = 1:iterations % for higher number of iteration, the result is very similar
   [q,r] = qr(A);
   A = q'*A*q; 
end

A, diag(A) % display results

I expected eigenvalues on diagonal of matrix A. Other numbers (outside of diagonal) should be 0.

Some of eigenvalues were correct (111 (dominant) and 0.0 ), but some are not correct.

The result is same for higher number of iteration (iterations = 80000).

It looks, that QR method works for dominant or assymetric eigenvalues ??

(I tested same QR algorithm for original matrix M, but with same incorrect result).

How can I fix it, for getting complete correct solutions ? Could you suggest me some suitable method / algorithm ?

Thank you

[Torsten]

https://www.andreinc.net/2021/01/25/computing-eigenvalues-and-eigenvectors-using-qr-decomposition

https://en.wikipedia.org/wiki/Eigenvalue_algorithm

Answer 1

[Nelson Rufus]

Hello! Your implementation of the QR algorithm for computing eigen values iteratively seems correct.

The matrices M, hM and A are singular (i.e. they have a zero eigenvalue). To see if that's affecting the QR algorithm, I tried shifting the M before it passes through the QR iterations and unshifting when results are displayed in the last line

M = magic(6);    % some sample matrix
shift = 0.2 .* eye(6); % the 0.2 is arbitrary
ideal = eigs(M)  % get ideal solution by original matlab function
M = M + shift;
% the implementation of QR
% ...

(A-shift), diag(A-shift) % display results

I found that this shifting does indeed resolve the discrepancy you were observing. The open question here is that when I shift using the following matrix,

shift = -111 .* eye(6);  % 111 is an eigen value of M!!

the QR algorithm still worked in spite of (M+shift) being singular.

A second approach to resolve your problem is to use shifted QR. This introduces a shift inside the QR loop. A crude implementation of this is as follows:

% using QR algorithm for detection eigenvalues
for k = 1:iterations % for higher number of iteration, the result is very similar
   local_shift_coeff = 0.0;
   if abs(A(3,2))>1.0e-8
       local_shift_coeff = A(3,2)
   end
   local_shift = local_shift_coeff * eye(6);
   [q,r] = qr(A-local_shift);
   A = r*q + local_shift;   
end

Using this approach, the QR algorithm solution does indeed match the eigenvalues using eigs(M). Notice above, that I'm checking only A(3,2) and accordingly picking the shift (local_shift_coeff). This was motivated by my observation that A(3,2) in your implementation was not decaying with iterations. A more general implementation should check all A(p+1,p) entries. (See https://www.math.usm.edu/lambers/mat610/class0331.pdf for an explanation of Shifted QR).