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?

168 次查看（过去 30 天）

ma Jack 2024-12-30，8:58

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2172557-hello-everyone-i-have-a-non-hermitian-matrix-but-i-find-that-even-using-vpa-i-cannot-obtain-the

评论： ma Jack about 18 hours 前

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.

21 个评论
显示 19更早的评论隐藏 19更早的评论

ma Jack 2024-12-31，13:47

在 MATLAB Online 中打开

TFig.fig

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 2024-12-31，21:24

编辑：Torsten 2025-1-1，2:56

在 MATLAB Online 中打开

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 2024-12-31，22:53

编辑：Paul 2025-1-1，14:44

在 MATLAB Online 中打开

eigtest.mat

@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 2025-1-4，7:37

在 MATLAB Online 中打开

VAPtest.fig

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