accuracy of eigenvalue calculation

I have this problem of finding eigenvector of some 6 by 6 matrix with the known 6 eigenvalues (lambda1, lambda2, etc). The form of the matrix is given as A=A0+D, where A0 is completely known and D is a diagonal matrix whose elements are not known.

The way I solved this was first I find the numerical solution to 6 characteristic eqautions of A with 6 diagonal elements (of D in A) as unknowns, i.e.,

det(A-lambda1*1)=0, det(A-lambda2*1)=0, etc.

Using fsolve, the numbers for diagonals were obtained.

Then I obtained the eigenvectors and eigenvalues of A, using eig function. While both eigenvectors and eigenvalues returned, eigenvalues obtained in this way is not the same as lambda1, lambda2, etc I started with. So I can not trust eigenvectors either. Why the eigenvalues are different? Shoudn't they be the same?

[Christine Tobler]

I don't think I understand what you're doing exactly. Is the matrix A you pass to eig symbolic or plain numeric with some choices for the diagonal values of D? Could you attach your code?

What does it mean that the eigenvalues of A=A0+D are known, but the values of D aren't known? Is the goal to find diagonal values of D that will lead to specific eigenvalues lambda1, ...?

[Gunn]

The goal of the first step is to find the diagonal values of D such that the eigenvalues of A=A0+D are the given lambda's. Here A0 and lambda's are all known.

Then as a second step, we pass A, whose elements are all known now, to a numeric eig function to obtain its eigenvalues and eigenvectors. Actually, I need to know the eigenvectors only. But if I want to make sure the eigenvectors are correct, as a consistency check, the eigenvalues obtained by eig function must be the same as lambda's. They are not and I do not undertsand why.

A0=[ 0.9996   -0.0078   -0.0000    0.0001    0.0000   -0.0000;
   -0.0074    0.9923   -0.0076    0.0000   -0.0001   -0.0001;
    0.0001   -0.0076    0.9922   -0.0075    0.0000   -0.0000;
    0.0001    0.0000   -0.0076    0.9923   -0.0076   -0.0000;
    0.0001   -0.0000   -0.0001   -0.0074    0.9923   -0.0076;
    0.0001    0.0000   -0.0001    0.0000   -0.0073    1.0000];
lambda_def=[0.9799;0.9857;0.9934;1.0011;1.0068;1.0089];
  
Id=eye(N);
  F=@(e)[det(A0+diag(e)-lambda_def(1)*Id); det(A0+diag(e)-lambda_def(2)*Id); det(A0+diag(e)-lambda_def(3)*Id); det(A0+diag(e)-lambda_def(4)*Id); det(A0+diag(e)-lambda_def(5)*Id); det(A0+diag(e)-lambda_def(6)*Id)];
  e0=[0.0001 ; 0.000001 ;-0.0007 ;-0.0003; 0.00008; 0.00005];
  [e1,fval,info]=fsolve(F,e0) 
  [V1, D1]=eig(A0+diag(e1))
  diag(D1) % must be the same as lambda_def

Answer 1

[Bruno Luong]

"Why the eigenvalues are different? Shoudn't they be the same? "

Because solving for e of equation using determinant

det(A0+diag(e)-lambda_def(i)*Id = 0, i=1,...,6

is a very bad numerical method to enforce A0+diag(e) to have eigen values of lambda_def, contrary to what they teach you in school.

I let you observe if lambda_def or lambda (obtained from eig) are better from this modified code (it make the same calculation, I just redefine F so I can use for different lambda):

A0=[ 0.9996   -0.0078   -0.0000    0.0001    0.0000   -0.0000;
   -0.0074    0.9923   -0.0076    0.0000   -0.0001   -0.0001;
    0.0001   -0.0076    0.9922   -0.0075    0.0000   -0.0000;
    0.0001    0.0000   -0.0076    0.9923   -0.0076   -0.0000;
    0.0001   -0.0000   -0.0001   -0.0074    0.9923   -0.0076;
    0.0001    0.0000   -0.0001    0.0000   -0.0073    1.0000];
lambda_def=[0.9799;0.9857;0.9934;1.0011;1.0068;1.0089];

N=size(A0);
Id=eye(N);
F=@(e,lambda) arrayfun(@(lambda) det(A0+diag(e)-lambda*Id),lambda);
e0=[0.0001 ; 0.000001 ;-0.0007 ;-0.0003; 0.00008; 0.00005];
[e1,fval,info]=fsolve(@(e) F(e,lambda_def),e0);
[V1, D1]=eig(A0+diag(e1));
lambda = diag(D1);

% Compute det criteria for input lambda_def and eigen-value from MATLAB EIG
F(e1,lambda_def) % <= LOOK CAREFUL THIS
F(e1,lambda)     % <= COMPARED TO THIS

Ask yourself which one is better, i.e. more accurate, eigen values of A+diag(e1).

[Gunn]

Thanks, Bruno.

But the part of the original problem I had was to find a diagonal perturbation to A0 with the initially given data of eigenvalues (lambda_def) and unperturbed matrix A0. Then the eventual goal was to find the eigenvectors of A0+D, which I guess can be easily obtained via eig function. But to confirm that the eigenvector obtained from A0+D via eig function is indeed a correct one, I needed to check on the eigenvalues by comparing those by eig function and the given values (lambda_def) I started with.

So I really need to find an accurate set of diagonals by solving the equations. Any suggestions to get these diagonals from the initially given data only?

[Bruno Luong]

Then use EIG, a professional code that the whole humanity uses since 60 years, instead your det formula.

A0=[ 0.9996   -0.0078   -0.0000    0.0001    0.0000   -0.0000;
   -0.0074    0.9923   -0.0076    0.0000   -0.0001   -0.0001;
    0.0001   -0.0076    0.9922   -0.0075    0.0000   -0.0000;
    0.0001    0.0000   -0.0076    0.9923   -0.0076   -0.0000;
    0.0001   -0.0000   -0.0001   -0.0074    0.9923   -0.0076;
    0.0001    0.0000   -0.0001    0.0000   -0.0073    1.0000];
lambda_target=[0.9799;0.9857;0.9934;1.0011;1.0068;1.0089];

opts = optimoptions(@fmincon, ...
    'StepTolerance', 1e-12, ...
    'OptimalityTolerance', 1e-16);
e0 = zeros(length(A0),1);
e = fmincon(@(e) lsqeig(e, A0, lambda_target), e0, [],[],[],[],[],[],[],opts);  

% Check how better eigen value match
eig(A0 + diag(e))
dlambda0 = deig(e0, A0, lambda_target)
dlambda = deig(e, A0, lambda_target)

%%
function dlambda = deig(e, A0, lambda_target)
lambda = eig(A0 + diag(e));
dlambda = sort(lambda)-sort(lambda_target);
end

%%
function f = lsqeig(e, A0, lambda_target)
dlambda = deig(e, A0, lambda_target);
f = sum(dlambda.^2);
end