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?

Question

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

评论： Torsten about 2 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.

11 个评论
显示 9更早的评论隐藏 9更早的评论

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

Paul 2025-1-1，14:55

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

在 MATLAB Online 中打开

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 about 2 hours 前

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.

请先登录，再进行评论。

请先登录，再回答此问题。

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?

11 个评论
显示 9更早的评论隐藏 9更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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?

11 个评论 显示 9更早的评论隐藏 9更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

11 个评论
显示 9更早的评论隐藏 9更早的评论