Why are the components of the two eigenvectors not continuous and smooth, when the underlying data is continuous and smooth?

For my research, I measure a 2x2 scattering (S) matrix for a number of frequency points. Each component of the S matrix is complex, but varies smoothly and continuously. I am interested in the inner product of the two eigenvectors of the S matrix, but it is discontinuous due to the components of the eigenvectors behaving oddly. The issue isn't that the label of which eigenvector is which has swapped, or a phase ambiguity of the eigenvectors. I define the eigenvectors as R1 = (A;B), R2 = (C;D) and choose A and D to be purely real. When the discontinuities in the inner product occur, a few things simultaneously occur. These include the label of eigenvalue 1 and 2 switching, a linear kink in the minima of A and D, a kink and switching of the real parts of B and C, and a kink or discontinuity in the imaginary parts of B and C. The code below and attached data show this behavior. Is there any way to fix this such that both the real and imaginary parts of the eigenvectors stay continuous, and the eigenvalues also be continuous with no switching?

%Experimental Data
load('Representative_Data.mat')

%Eigensolver
for k=1:length(x_pts)
    [Right_evec(:,:,k),eval(:,:,k)] = eig(smat(:,:,k));
    Eval_1(k) = eval(1,1,k);
    Eval_2(k) = eval(2,2,k);
end

%Get rid of eigenvector phase ambiguities
Phase_Component_1 = angle(Right_evec(1,1,:)); %A
Phase_Component_2 = angle(Right_evec(2,2,:)); %D
Right_evec(:,1,:) = Right_evec(:,1,:).*exp(-1i*Phase_Component_1);
Right_evec(:,2,:) = Right_evec(:,2,:).*exp(-1i*Phase_Component_2);

for k=1:length(x_pts)
    Inner_Product(k) = Right_evec(:,1,k)'*Right_evec(:,2,k);
end

Marker_size = 20;

%Eigenvalues
figure('Position',[200,100,1500,730]);
scatter(x_pts,real(Eval_1),Marker_size,'filled'); hold on;
scatter(x_pts,imag(Eval_1),Marker_size,'filled');
scatter(x_pts,real(Eval_2),Marker_size,'filled');
scatter(x_pts,imag(Eval_2),Marker_size,'filled');
xlabel('Frequency (GHz)')
ylabel('Eigenvalues')
set(gca,'Fontsize',26)
legend('Re \lambda_1','Im \lambda_1','Re \lambda_2','Im \lambda_2')
grid on; box on;

%Eigenvector Components A and D
figure('Position',[200,100,1500,730]);
scatter(x_pts,real(squeeze(Right_evec(1,1,:))),Marker_size,'filled'); hold on;
scatter(x_pts,real(squeeze(Right_evec(2,2,:))),Marker_size,'filled');
%scatter(x_pts,imag(squeeze(Right_evec(1,1,:))),Marker_size,'filled'); %Imaginary part specified to be 0
%scatter(x_pts,imag(squeeze(Right_evec(2,2,:))),Marker_size,'filled'); %Imaginary part specified to be 0
xlabel('Frequency (GHz)')
ylabel('Eigenvector components')
set(gca,'Fontsize',26)
legend('Re A','Re D')
grid on; box on;

%Eigenvector Components B and C
figure('Position',[200,100,1500,730]);
scatter(x_pts,real(squeeze(Right_evec(2,1,:))),Marker_size,'filled'); hold on;
scatter(x_pts,imag(squeeze(Right_evec(2,1,:))),Marker_size,'filled');
scatter(x_pts,real(squeeze(Right_evec(1,2,:))),Marker_size,'filled');
scatter(x_pts,imag(squeeze(Right_evec(1,2,:))),Marker_size,'filled');
xlabel('Frequency (GHz)')
ylabel('Eigenvector components')
set(gca,'Fontsize',26)
legend('Re B','Im B','Re C','Im C')
grid on; box on;

%Inner Product
figure('Position',[200,100,1500,730]);
scatter(x_pts,real(Inner_Product),Marker_size,'filled'); hold on;
%scatter(x_pts,imag(Inner_Product),Marker_size,'filled'); %Imaginary part is properly continuous
xlabel('Frequency (GHz)')
ylabel('Re <R_1|R_2>')
set(gca,'Fontsize',26)
grid on; box on;

[David Goodmanson]

Hi Jared,

With eig, there is no guarantee about the order of eigenvalues and their associated eigenvectors. Each eigenvector does have norm 1 and that's all that is guaranteed.

It looks like eig has long stretches that preserve continuity and it occasionally swaps the order of eigenvalues from the order that would be implied by continuity. So you could find discontinuities in say, the (1,1) eigenvalue as a function of frequency and at that point swap the eigenvalues given by eig (and swap the associated eigenvector columns at the same time), keep them swapped until you find the next discontinuity and then swap them back, etc.

WIth eigenvectors you can multiply by a phase factor which you are doing by setting the first component to be real. However, first-component-real is an arbitrary choice unless there is some physics/ EE reason to do so, and those phase factors persist into the calculation of Re(<R1|R2>). So it's not clear what the conclusion will be. This is in contrast to the phases contained in the S matrix itself, which are significant.

0 Comments

[Jared Erb]

Hi David,

Thanks for responding. I understand that there's no guarantee about the order, and that is ok and easy to fix just as you mentioned. For the inner product though the issue isn't just that the order is switching, if that were the case it would just get a strict negative sign at the discontinuities. Making the first components real is arbitrary and unnecessary, but it makes the discontinuities a bit less extreme. Since the discontinuities in the inner product aren't entirely due to the order switching, I was wondering if there was a different way to fix the discontinuities, at least up to a minus sign. Because even if the order is fixed by changing the eigenvalues to be continuous, the inner product would still be discontinuous. It would also be interesting to know why the order switches when it does, but that might not be easy to determine. Thank you very much.

[Jared Erb]

If I choose the eigenvector components A and C to be purely real (instead of A and D), and fix the eigenvalues to be continuous (fixing the eigenvectors at the same time), then the inner product (and eigenvector components) becomes properly continuous.

Answer 1

[Matt J]

Note that eigenvectors are not generally continuous functions of the matrix entries, as shown in this old example from J.W. Givens:

H(0)=eye(2)

H(k)=[1 + k*cos(2/k),  -k*sin(2/k);
     -k*sin(2/k),      1 - k*cos(2/k)]

The matrix H(k) is continuous as a function of k, even at k=0, but its matrix of eigenvectors U is not:

 U(k)=[sin(1/k) -cos(1/k);
       cos(1/k)  sin(1/k)]

I believe this can only occur in a neighborhood of H where eigenvalues have multiplicity greater than 1. Based on your plots, it is not clear whether this happens for your data. In general, though, it should be kept in mind that the eigenvalues of a matrix, but not necessarily the eigenvectors, are continuous in the matrix elements.

Answer 2

[John D'Errico]

You can think of the eigenvalues (and assiciated eigenvectors) as coming from the roots of a polynomial. In fact, they do exactly that. They are derived from a characteristic polynomial. Now, suppose you have a polynomial where the coefficients are a function themselves of some parameter. The roots will vary, as a function of that parameter. The problem is, those roots (as well as the function roots) do not understand the parameter they depend upon. That are just numbers. The function roots does not understand that parameter. It gives not a tinker's damn why it is called, only that you pass it a vector of coefficients, and it returns a valid result.

For example, suppose we consider a simple polynomial in x. I'll just pick something at random. Again, while this is a polynomial, I could as easily be creating a matrix, and computing the eigenvalues. There is a one to one correspondence here. We can always create a matrix that has the same eigenvalues as any polynomial (this is called the companion matrix, created by the matlab function compan), and we can always take any matrix and generate the characteristic polynomial. So what I do here is entirely valid as a comparison. I'll show it in symbolic form first.

syms x t
P(x) = x^5 - t/2*x^4 + 3*x^3 + (t^2 + 1)*x^2 + 2*t*x - t^3

P(x) = 

I've been very arbitrary in that polynomial, and honestly, I have not a clue what it looks like, as a function of t, and where the roots might lie, nor how they might vary.

coefflist = coeffs(P,x)

coefflist(x) = 

Hmm. Coeffs decided to put them in the opposite order of what roots would want., and for some stupid reason, flip does not fix the issue. I care not, and I won't bother fixing it. ;-)

Pcoeffs = @(t) [1, -t/2 ,3 , t^2+1 ,2*t, -t^3]

Pcoeffs = function_handle with value:
    @(t)[1,-t/2,3,t^2+1,2*t,-t^3]

For any value of t, we can comput the roots, and follow them, as they are returned from roots. Again, roots does not care about the order of those roots, or how that order is related to the parameter t. You NEED to understand that. It does not think of the roots as a function of t. Only you think that way. and roots gives not a hoot about what you understand.

nr = 1001;
tlist = linspace(-10,10,nr);
rlist= zeros(nr,5);
for i = 1:nr
  rlist(i,:) = roots(Pcoeffs(tlist(i))).';
end

Now, plot the roots, in sequence, exactly as they were returned.

plot(real(rlist),imag(rlist),'-')

Do you see the roots bounce around in order? If, instead, I resequence those roots, they would surely follow nice paths in the complex plane. Using my eigenshuffle code (Found on the file exchange, it can improve on that mess, though even it is not perfect.)

Aseq = zeros(5,5,nr);
for i = 1:nr
  Aseq(:,:,i) = compan(Pcoeffs(tlist(i)));
end
[Vseq,Dseq] = eigenshuffle(Aseq);
plot(plot(real(Dseq.'),imag(Dseq.'),'-'))

Here you can see the roots do indeed follow well-behaved paths around the complex plane, though eigenshuffle was not perfect in its sequencing choice.

Again, you need to understand that eig does not know how the eigenvectors and eigenvalues should vary as the elements of a matrix vary.