Hello everyone, I have a non-Hermitian matrix, but I find that even using vpa, I cannot obtain the correct eigenvalues. Does anyone know how to resolve this?

Below are my code and .mat data. The green circles represent the correct eigenvalue range for the non-Hermitian matrix (derived using a complex mathematical method; since the method is too complicated, I’m directly using the results obtained from it). The black points represent the eigenvalues computed using the vpa function, but they do not coincide with the green circles. (Indeed, the number of green circles and black points is different, but if the black points were correct, their upper and lower bounds should overlap.) This means that I didn’t obtain the correct eigenvalues. I know that non-Hermitian matrices cannot be exactly diagonalized, but I’ve heard from others that using the high-precision vpa function can help, but it doesn’t seem to work in my case.

load list_beta0.mat
load GBZband.mat

N=90;
t1=0.5;
t2=0.5;

t3=1/5;
gamma1=5/3;
gamma2=1/3;

HH=Hobc(t1,t2,t3,gamma1,gamma2,N);
[~,val]=eig(vpa(HH,34));
tval=diag(val);
[~,index]=sort(real(tval));

sort_valR=tval(index);

figure
hold on
scatter(t1.*ones(size(GBZband)), abs(GBZband),10,'go')
plot(t1.*ones(size(sort_valR)),abs(sort_valR),'k.','MarkerSize',8)
ggg=legend('GBZ-eigenval','OBC-eigenval');
set(ggg,'box','off')
hold off




function [op]=Hobc(t1,t2,t3,gamma1,gamma2,N)
a=[1;0];
b=[0;1];
H=zeros(2*N,2*N);
n1=zeros(N,1);
n2=zeros(N,1);

for n=1:N
    n1(n)=1;
    n2(n+1)=1;

    H=H+(t1+gamma1/2).*kron(n1,a)*kron(n1',b')+(t1-gamma1/2).*kron(n1,b)*kron(n1',a');

    if n~=N
        H=H+(t2+gamma2/2).*kron(n1,b)*kron(n2',a')+(t2-gamma2/2).*kron(n2,a)*kron(n1',b')+...
            t3.*(kron(n1,a)*kron(n2',b')+kron(n2,b)*kron(n1',a'));
    end

    n1=zeros(N,1);
    n2=zeros(N,1);

end

op=H;
end

Does anyone have any insights on this issue? I look forward to any responses. If the question is unclear, please let me know.

[Torsten]

The black and the green points share a common range (from 0 to 0.9, I guess) - and most probably there is a green point below each black point. So I don't understand what the problem is.

[ma Jack]

Thank you, sir, for your comments. First, I’d like to clarify that the green are not points, but circles, so they do not obscure the black points. In other words, the two markers do not overlap. Also, perhaps I provided too few examples. You can see from the following plot, which is drawn for a wide range of parameters, that there are both good and bad cases. For instance, in the blue box, the black points and green circles' upper and lower bounds match very well; whereas in the red box, while the upper bounds match very well, the lower bounds do not. My goal is to have both the upper and lower bounds match as closely as they do in the blue box, and this is the issue I am facing.

The mathematical definition of a Non-Hermitian matrix is H \neq H^{\dagger}, where H^{\dagger} is the conjugate transpose of the matrix H.

[Torsten]

The number of green circles is much larger than the number of black points. So maybe the lower bound would have been met better if you had chosen a higher dimension for HH.

But since we do know nothing about the problem you try to solve here, we can only speculate. I don't think that from the information given, you can expect a useful answer.

[Paul]

Why is vpa necessary? I ran the code and compared the vpa eig to the double eig and the eigenvalues laid on top of each other. Does the matrix HH possess some attribute(s) that renders Matlab's eig solutions ineffective?

[ma Jack]

Thank you all for your comments. It seems that my question wasn’t very clear. Let’s look at an equivalent problem instead. Below, I’ve provided a piece of code that uses two methods to compute the eigenvalues of the same matrix. I’ve plotted them on the same graph: black points are from the numerical method, and blue circles are from the analytical method (which is undoubtedly the correct result, but very time-consuming). As you can see, the part enclosed by the red box shows that the numerical solution and the analytical solution do not perfectly match (the mismatch becomes more severe as N increases). This is the issue I’m facing. I want to understand the reason for this mismatch and how to adjust the numerical computation to make it disappear. Below are my code and the .fig format image.

clear;clc;close all

w=1;
u=2.5;

N=30;


valsto=zeros(2*N+2,100);
coun=0;
for v1=linspace(-3,3,100)
    coun=coun+1;

    [~,val]=eig(H(u,v1,w,N));
    tval=diag(val);
    [~,index]=sort(real(tval));
    val=tval(index);
    valsto(:,coun)=val;

end

figure
hold on
plot(linspace(-3,3,100),abs(valsto(:,:)),'k.')

%%%%%%%%%%=========%%%%%%%%%%%%%%
%%%%%%%%%%=========%%%%%%%%%%%%%%



clearvars -except N
syms u v w real


val2=eig(H(u,v,w,N));
val2=subs(val2,{u,w},{2.5,1});

valsto2=zeros(2*N+2,100);
coun=0;
for v1=linspace(-3,3,100)
    Mark=coun/100

    coun=coun+1;


    tval=val2;
    tval=double(subs(tval,v,v1));
    [~,index]=sort(real(tval));
    tval=tval(index);
    valsto2(:,coun)=tval;

end

plot(linspace(-3,3,100),abs(valsto2(:,:)),'bo')
hold off


function [H0]=H(u,v,w,N)

T1=[];
T2=[];
for hh=1:N
    T1=[T1;v+u;w+u];
    T2=[T2;v-u;w-u];
end
T1=[T1;v+u];
T2=[T2;v-u];
H0=diag(T1,1)+diag(T2,-1);


end

[Torsten]

It seems that the numerical "eig" is not able to compute the (purely imaginary) eigenvalues of the matrix correctly.

I have no experience with eigenvalue problems - so I cannot tell you the reason. Maybe it's because the eigenvalues are very close to each other.

Maybe @Bruno Luong has an idea.

But the ranges of the absolute values of the eigenvalues computed symbolically and numerically differ from your figure above at v = -2.5... - they look almost identical. What MATLAB version do you use ?

format long
vnum = linspace(-3,3,100);
N = 30;
unum = 2.5;
wnum = 1;
Lnum = eig(H(unum,vnum(10),wnum,N))

Lnum = 
 -0.000000000000004 + 0.000000637176401i
 -0.000000000000004 - 0.000000637176401i
  0.117913226667825 + 2.815209265412365i
  0.117913226667825 - 2.815209265412365i
  0.058263159695574 + 2.847351613723787i
  0.058263159695574 - 2.847351613723787i
 -0.003527084064689 + 2.857496369370038i
 -0.003527084064689 - 2.857496369370038i
 -0.065451748175237 + 2.845500743757776i
 -0.065451748175237 - 2.845500743757776i
 -0.125551134748330 + 2.811591954120779i
 -0.125551134748330 - 2.811591954120779i
  0.173282561441625 + 2.761452452122647i
  0.173282561441625 - 2.761452452122647i
 -0.181781668336172 + 2.756321445319460i
 -0.181781668336172 - 2.756321445319460i
  0.222183382929459 + 2.686234285694165i
  0.222183382929459 - 2.686234285694165i
 -0.232082229146546 + 2.680545724956049i
 -0.232082229146546 - 2.680545724956049i
  0.263746680442496 + 2.588555362681767i
  0.263746680442496 - 2.588555362681767i
 -0.274862100741577 + 2.585453886527108i
 -0.274862100741577 - 2.585453886527108i
  0.302890771046344 + 2.471955990603175i
  0.302890771046344 - 2.471955990603175i
  0.145849001350266 + 2.510232327098566i
  0.145849001350266 - 2.510232327098566i
 -0.309905715522997 + 2.473283236090422i
 -0.309905715522997 - 2.473283236090422i
  0.334249947896254 + 2.348418957711169i
  0.334249947896254 - 2.348418957711169i
 -0.087485083148853 + 2.430947123540911i
 -0.087485083148853 - 2.430947123540911i
 -0.336706094025641 + 2.349493259761486i
 -0.336706094025641 - 2.349493259761486i
  0.350119758606805 + 2.219641791931203i
  0.350119758606805 - 2.219641791931203i
 -0.350876110431463 + 2.219949033114284i
 -0.350876110431463 - 2.219949033114284i
  0.347483747230019 + 2.088634094777779i
  0.347483747230019 - 2.088634094777779i
 -0.347693883780828 + 2.088684401654675i
 -0.347693883780828 - 2.088684401654675i
  0.323575519219198 + 1.960727392120438i
  0.323575519219198 - 1.960727392120438i
 -0.323624238534635 + 1.960730968148195i
 -0.323624238534635 - 1.960730968148195i
  0.275937695315164 + 1.843213722793182i
  0.275937695315164 - 1.843213722793182i
 -0.275946931262473 + 1.843214126147148i
 -0.275946931262473 - 1.843214126147148i
 -0.203496648292608 + 1.745683177169336i
 -0.203496648292608 - 1.745683177169336i
  0.203495211261072 + 1.745682622497439i
  0.203495211261072 - 1.745682622497439i
 -0.108701766122748 + 1.679727642760095i
 -0.108701766122748 - 1.679727642760095i
  0.108701703664981 + 1.679727352055979i
  0.108701703664981 - 1.679727352055979i
  0.000000069567721 + 1.656185300634534i
  0.000000069567721 - 1.656185300634534i

%%%%%%%%%%=========%%%%%%%%%%%%%%
syms u v w real

H(u,v,w,N)

ans = 

Lsym=vpa(eig(subs(H(u,v,w,N),[u v w],[unum,vnum(10),wnum])))

Lsym = 

function [H0]=H(u,v,w,N)

T1=[];
T2=[];
for hh=1:N
    T1=[T1;v+u;w+u];
    T2=[T2;v-u;w-u];
end
T1=[T1;v+u];
T2=[T2;v-u];
H0=diag(T1,1)+diag(T2,-1);


end

[Paul]

@Torsten

Interesting that the problem has changed from the original question. The H matrix is different and now we are comparing double and sym eig, whereas the original question was comparing double eig to a "complex mathematical method."

As formulated now, for the second case the eigenvalues should be computed symbolically first in terms of u,v,w and then evaluated for the specific values of u,v1, and w.

%Lsym=vpa(eig(subs(H(u,v,w,N),[u v w],[unum,vnum(10),wnum])))
% should be
%Lsym=vpa(subs(eig(H(u,v,w,N)),[u v w],[unum,vnum(10),wnum]))

When I did that using vnum(9) (not vnum(10) as in your example), the eigenvalues have real parts as well and they are different than using double from the outset.

load eigtest eigsym v1 HHs HHd
HHs(1:5,1:5)

ans = 

unum = 2.5;
wnum = 1;
norm(HHd - double(subs(HHs,[sym('u'),sym('v'),sym('w')],[unum,v1,wnum])))

ans = 1.3878e-17

doubleeig = eig(HHd);
symeig    = subs(eigsym,[sym('u'),sym('v'),sym('w')],[unum,v1,wnum]);
dsymeig = double(vpa(symeig));
figure
plot(real(doubleeig),imag(doubleeig),'o',real(dsymeig),imag(dsymeig),'x')
axis equal

Note that the double solution yields two eigenvalues very close to the origin, whereas the sym->vpa->double solution only yields one.

The symbolic eigenvalues are expressed as square roots of roots of a 31st order polynomial. I wonder about the accuracy of the solutions of that polynomial after converting to double as above. The following plot is an attempt to compare the double solution to the sym->vpa->double solution.

figure
hold on
for ii = 1:numel(doubleeig)
    plot(ii,min(svd(HHd-doubleeig(ii)*eye(size(HHd)))),'bo');
end
for ii = 1:numel(dsymeig)
    plot(ii,min(svd(HHd-dsymeig(ii)*eye(size(HHd)))),'rx');
end

[ma Jack]

@Torsten It doesn't seem to be an issue with the MATLAB version, as I am using MATLAB R2023b. This problem likely exists more generally.

[Torsten]

I have no experience with eigenvalue problems - so I cannot tell you the reason. Maybe it's because the eigenvalues are very close to each other.

Sorry.

[Paul]

Consider a smaller problem, with N = 6

N = 6;
syms u v w s
H0 = H(u,v,w,N);

The characterisic polynomial of H0 is

c = collect(charpoly(H0,s),s);

Let v = -u and solve for the roots

v = -u;
e = solve(subs(c),s)

e = 

We see two repeated eigenvalues at the origin, and two sets of six repeated eigenvalues. For the prescribed values of u and w these eigenvalues are imaginary.

I'm going to speculate that it becomes very difficult to solve for the eigenvalues of H0 in double precision as N gets larger and v approaches -u (for the prescribed values of u and w), which I think is basically the problem of solving for near-repeated roots of a high order polynomial (note that roots actually solves for the roots using eig).

@Christine Tobler can probably provide more insight into solving for the eigenvalues of this particular matrix.

function [H0]=H(u,v,w,N)
T1=[];
T2=[];
for hh=1:N
    T1=[T1;v+u;w+u];
    T2=[T2;v-u;w-u];
end
T1=[T1;v+u];
T2=[T2;v-u];
H0=diag(T1,1)+diag(T2,-1);
end

[Torsten]

I agree with @Paul. At v = -2.5, your matrix has two N-fold eigenvalues. In the neighborhood of v = -2.5, these eigenvalues spread to N distinct eigenvalues each, but it will be very difficult for the solver to separate them because they lie very close to each other.

[Paul]

Using the generalized eigenvalue solver, as previously discussed here, doesn't seem to make much of a difference, but it does seem to get rid of that lone pair of double eigenvalues that are very close to the sym solution.

load eigtest eigsym v1 HHs HHd
unum = 2.5;
wnum = 1;
doubleeig1 = eig(HHd);
doubleeig2 = eig(HHd,eye(size(HHd)),'chol');
doubleeig3 = eig(HHd,eye(size(HHd)),'qz');

symeig    = subs(eigsym,[sym('u'),sym('v'),sym('w')],[unum,v1,wnum]);
dsymeig = double(vpa(symeig));
figure
hold on
plot(real(doubleeig1),imag(doubleeig1),'o','DisplayName','deig');
plot(real(dsymeig),imag(dsymeig),'x','DisplayName','symeig');
axis equal
hold on
plot(real(doubleeig2),imag(doubleeig2),'d','DisplayName','geig-chol');
plot(real(doubleeig3),imag(doubleeig3),'s','DisplayName','geig-qz');
legend

●

[ma Jack]

Thank you for all the comments on my question. I believe I’ve solved the issue—actually, it’s quite simple and doesn’t require complex discussions. The solution is to use the vpa function, and it’s universal. Below is the code:

clear;clc;close all
tic

delete(gcp('nocreate'));
c=parcluster('local');
c.NumWorkers=50;
parpool(c, c.NumWorkers);


w=1;
u=2.5;

N=30;


valsto=zeros(2*N+2,100);
v10=linspace(-3,3,100).';


parfor jj=1:length(v10)


    v1=v10(jj);

    [~,val]=eig(vpa(H(u,v1,w,N),30));
    %[~,val]=eig(H(u,v1,w,N));
    tval=diag(val);
    [~,index]=sort(real(tval));
    val=tval(index);
    valsto(:,jj)=val;

end

figure
hold on
plot(linspace(-3,3,100),abs(valsto(:,:)),'k.')

%%%%%%%%%%=========%%%%%%%%%%%%%%
%%%%%%%%%%=========%%%%%%%%%%%%%%



clearvars -except N
syms u v w real


val2=eig(H(u,v,w,N));
val2=subs(val2,{u,w},{2.5,1});

valsto2=zeros(2*N+2,100);
v10=linspace(-3,3,100).';


parfor jj=1:length(v10)


    v1=v10(jj);
    tval=val2;
    tval=double(subs(tval,v,v1));
    [~,index]=sort(real(tval));
    tval=tval(index);
    valsto2(:,jj)=tval;

end

plot(linspace(-3,3,100),abs(valsto2(:,:)),'bo')
hold off
%saveas(gcf,"VAPtest0.fig")
toc

function [H0]=H(u,v,w,N)

T1=[];
T2=[];
for hh=1:N
    T1=[T1;v+u;w+u];
    T2=[T2;v-u;w-u];
end
T1=[T1;v+u];
T2=[T2;v-u];
H0=diag(T1,1)+diag(T2,-1);


end

The attachment contains the plot of the eigenvalues obtained using vpa computation and symbolic computation.

[Torsten]

Is the method how MATLAB internally computes the eigenvalues different for the vpa and the symbolic computation ? I'm not sure.

[Paul]

I ran the udpated code in this comment and got the same eigenvalues for the vpa solution as shwon on the plot in this comment for jj = 9. I thought that throughout this thread the vpa eig solution was being taken as "truth" and the question was whether or not the double eig solution could be as accurate. So I'm not sure what issue has been solved (unless coverting to vpa then calling eig is a better solution than calling eig then converting to vpa?); can you clarify what goal has been achieved? Of course, if you're satisfied that's all that counts.

[Torsten]

@Paul

Maybe

[~,val]=eig(vpa(H(u,v1(9),w,N),30));

uses the "normal" numerical method to determine eigenvalues (only with variable-precision arithmetic) while

val2=eig(H(u,v,w,N));
val2=subs(val2,[u,v,w],[2.5,v1(9),1]);

directly tries to find the roots of the characteristic polynomial using "root" ?

At least, the first method seems to be much faster than the second.

That's why I wondered about the difference between the two methods used in this comment .

[Christine Tobler]

I don't know much about the symbolic version of eig. For numerics, yes we wouldn't compute eig through the characteristic polynomial, but when there are symbolic variables I don't think eig has an option. So it's very possible eig for vpa (variable-precision numbers, but no variables where values aren't known) is using a separate algorithm.

In terms of the numerics, if you call max(H(u,v1,w,N)) in the "numerics" loop above, you'll see some of these eigenvalues have condition numbers up to 1e40, which could easily explain the problems seen in that computation.

[Torsten]

@Christine Tobler

Is the "usual" numerical "eig" function able to handle variable-precision arithmetics input matrices as in the command

[~,val]=eig(vpa(H(u,v1(9),w,N),30));

and return solutions of the same precision ? Or is there a switch of the algorithm in this case ?

What about other numerical functions if inputs have enhanced precision ?

[Christine Tobler]

Sorry, that was unclear. The variable-precision case does not flow through the "usual" numerical eig function, all algorithms for vpa in Symbolic Toolbox are separate from the ones for numeric matrices in core MATLAB.

I know very little about Symbolic Toolbox, so this was all just speculation. From a numerics perspective, while we can convert eigenvalue problems to polynomial problems and vice-versa, both will be different in terms of error analysis.

There, we use eig even to find roots of a polynomial (see roots.m); but I don't think these perspectives translate to the symbolic or variable-precision case.

[Torsten]

@Christine Tobler

Thank you.