How to avoid using det, when looking for the complex root w with det(M(w)) = 0

Hi all,

I want to find a complex number w such that det(M(w)) converges to 0 (note that M is a matrix and it is a function of w), but the det function seems to have a large error, how can I avoid using it?M is an 8*8 matrix, so it would be very complicated to write out its determinant expression.

Thanks in advance.

[KSSV]

Show us your M...

[ma Jack]

GvXY is the matrix M,Zero_seek_Main() is the main function

function []=Zero_seek_Main()
clear;clc;close all

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

d=500*1e-9;
rx=d/4;
ry=rx;
num_Re=5;
num_Im=5;
%====================

%==========%
lsqoptions = optimoptions(@lsqnonlin,'Display','none');
%==========%

Real_w0=linspace(0.8,1.2,num_Re);
Imag_w0=linspace(-0.2,0.2,num_Im);


Result_stoXY=zeros(num_Re*num_Im,1);
count=0;

for uu=1:num_Re
        
        for vv=1:num_Im
            count=count+1;
            
            mark=count/(num_Re*num_Im)
            
            Re_w=Real_w0(uu);
            Im_w=Imag_w0(vv);
            w=Re_w+1j*Im_w;
            
            
            FUNXY=@(x)det(GvXY(x,d,rx,ry));

            
            %======================
            [find_x,~,~,~]=lsqnonlin(FUNXY,w,[],[],lsqoptions);
            Result_stoXY(count,1)=find_x;
            %======================
        end
end

     
figure
Result_index=[1:length(Result_stoXY)].';
plot(Result_index,Result_stoXY,'bo')
            
            
            
end

function [ModeXY]=GvXY(w,d,rx,ry)
c0=299792458;%light
wsp=3.742*1e15;%rad/s
eps0=8.8541*1e-12;


k=w*wsp/c0;
M1=RegSigGtXY(k,d);%size:2*2
M2=SigGtXY(k,rx,ry,d);%size:2*2


ModeXY=[M1,M2,M2,M2;...
                  M2,M1,M2,M2;...
                  M2,M2,M1,M2;...
                  M2,M2,M2,M1];%size:(2*4,2*4)

ModeXY=(k^2)/eps0.*ModeXY;



ModeXY(isnan(ModeXY))=0;%set NAN to 0
ModeXY(isinf(ModeXY))=0;%set INF to 0

%ModeXY=ModeXY(1:4,1:4);%If I enable this command. Then lsqnonlin will not 
%report an error, if I increase the 
%size of the matrix (for example, 5*5) lsqnonlin will report an error

end

function [o4]=RegSigGtXY(k,d)

N=1;

index_sto=zeros((2*N+1)^2-1,2);
cc=0;
for m=-N:N
    for n=-N:N
        
        if m==0&&n==0
         cc=cc;
        else
        cc=cc+1;
        
        index_sto(cc,1:2)=[m,n];
        end
        
    end
end

Sum_sto11=zeros((2*N+1)^2-1,1);%xx
Sum_sto22=zeros((2*N+1)^2-1,1);%yy
Sum_sto12=zeros((2*N+1)^2-1,1);%xy
Sum_sto21=zeros((2*N+1)^2-1,1);%yx

parfor ii=1:size(index_sto,1)
m=index_sto(ii,1);
n=index_sto(ii,2);

%==================
nR=abs(m*d+1j*n*d);
Rx=m*d;
Ry=n*d;
B=BB(k,nR);
A=AA(k,nR);

term=exp(1j*k*nR)/(4*pi*nR);
%==================


Sum_sto11(ii,1)=B/(nR^2)*(Rx^2)*term+A*term;
Sum_sto22(ii,1)=B/(nR^2)*(Ry^2)*term+A*term;
Sum_sto12(ii,1)=B/(nR^2)*Rx*Ry*term;
Sum_sto21(ii,1)=B/(nR^2)*Ry*Rx*term;
end

o4=zeros(2,2);

o4(1,1)=sum(Sum_sto11);
o4(2,2)=sum(Sum_sto22);
o4(1,2)=sum(Sum_sto12);
o4(2,1)=sum(Sum_sto21);


end

function [o3]=SigGtXY(k,rx,ry,d)

N=1;



index_sto=zeros((2*N+1)^2,2);
cc=0;
for m=-N:N
    for n=-N:N
        cc=cc+1;
        
        index_sto(cc,1:2)=[m,n];
    end
end

Sum_sto11=zeros((2*N+1)^2,1);%xx
Sum_sto22=zeros((2*N+1)^2,1);%yy
Sum_sto12=zeros((2*N+1)^2,1);%xy
Sum_sto21=zeros((2*N+1)^2,1);%yx

parfor ii=1:size(index_sto,1)
m=index_sto(ii,1);
n=index_sto(ii,2);

%=====================
nR=abs(m*d+n*d*1j+rx+ry*1j);
Rx=m*d+rx;
Ry=n*d+ry;
Term=exp(1j*k*nR)/(4*pi*nR);
A=AA(k,nR);
B=BB(k,nR);
%=====================

Sum_sto11(ii,1)=B/(nR^2)*(Rx^2)*Term+A*Term;
Sum_sto22(ii,1)=B/(nR^2)*(Ry^2)*Term+A*Term;
Sum_sto12(ii,1)=B/(nR^2)*Rx*Ry*Term;
Sum_sto21(ii,1)=B/(nR^2)*Ry*Rx*Term;
end

o3=zeros(2,2);

o3(1,1)=sum(Sum_sto11);
o3(2,2)=sum(Sum_sto22);
o3(1,2)=sum(Sum_sto12);
o3(2,1)=sum(Sum_sto21);

end

function [o1]=AA(k,R)
o1=1+(1j*k*R-1)/((k*R)^2);
end

function [o2]=BB(k,R)
o2=(3-3*1j*k*R-(k*R)^2)/((k*R)^2);
end

[KSSV]

What exactly you are trying to do? x is not defined in your case.

[ma Jack]

x is w, which is defined in the following code:

FUNXY=@(x)det(GvXY(x,d,rx,ry));

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

[find_x,~,~,~]=lsqnonlin(FUNXY,w,[],[],lsqoptions);

[Walter Roberson]

As discussed in your previous question, you have the difficulty that you are working with a matrix whose determinant is on the order of 10^250

On the first ModeXY matrix that is generated, there are only four unique values. If you construct a representative symbolic matrix in terms of the pattern of unique values, and take det() of it, and substitute in the unique values, then that is a lot faster than calculating det() of the original matrix symbolically. The idea of calculating it symbolically being to reduce the error in the calculation of det()

ModeXY = [M1,M2,M2,M2;...
                  M2,M1,M2,M2;...
                  M2,M2,M1,M2;...
                  M2,M2,M2,M1];%size:(2*4,2*4)

That promises that the pattern continues of there being only 4 unique elements in the matrix, so you can

pre-calculate the determinant as

(V1 + V2 - V3 - V4)^3*(V1 - V2 - V3 + V4)^3*(V1 + V2 + 3*V3 + 3*V4)*(V1 - V2 + 3*V3 - 3*V4)

which would be 0 if and only if any one of the sub-expressions is 0, which happens if

V1 + V2 = V3 + V4
V1 + V4 = V2 + V3
V1 + 3*V3 = V2 + 3*V4
V1 + V2 + 3*V3 + 3*V4 == 0

You might be able to take advantage of those to seek for a zero with a lower range.

[ma Jack]

Sir thank you for your suggestion, but I think Mr. Matt J's answer is better.

Answer 1

[Matt J]

Because your matrix appears to be symmetric, I suggest minimizing instead norm(M(w)) which is the maximum absolute eigenvalue of M. This is the same as forcing M to be singular.

EDIT: rcond(M) is probably more appropriate than norm(M)

Additionally, I suggest using fminsearch instead of lsqnonlin, since you only have a small number of variables and a non-differentiable cost function. Be mindful, however, that you must express your objective function in terms of a vector z of real variables.

w=@(z) complex(z(1),z(2));

zopt=fminsearch(@(z) norm(M(w(z)))  ,z0)

wopt=w(zopt)

[ma Jack]

Sorry, I don't quite understand why norm(M(w)) is used instead of det(M(w)), is this some kind of mathematical knowledge?

Also, is there any advantage of fminsearch over lsqnonlin? Why should we use it

[Matt J]

I don't quite understand why norm(M(w)) is used instead of det(M(w)),

John demonstrated the problems with det() in your earlier post. Let's see how norm(M) does in similar examples:

M1=eye(1200)/2;
M2=0*M1;

det(M1), det(M2)

ans = 0

ans = 0

norm(M1), norm(M2)

ans = 0.5000

ans = 0

Clearly norm(M) is better able to distinguish between singular and non-singular matrices than det. rcond() would be even better because it is independent of the scaling of M1:

rcond(M1),rcond(1e-100*M1), rcond(M2)

ans = 1

ans = 1

ans = 0

is there any advantage of fminsearch over lsqnonlin?

lsqnonlin assumes that the cost function is differentiable whereas fminsearch does not. It's not clear to me that your M(w) is a differentiable function of w, so fminsearch might therefore be safer. But you can also try lsqnonlin if you wish.

[Matt J]

ma Jack's comment moved here:

Thank you very much for your reply sir, but I think I need to spend some time to understand and test them.

[ma Jack]

Sorry sir, I have one more doubt, if we find a complex w such that norm(M(w)) -> 0,will this also make det(M(w)) -> 0? Please excuse my tardiness.

[Matt J]

Yes. All of the above are conditions for M --> singular matrix.

[Paul]

@Matt J

For a matrix M, norm(M) is the maximum singular value of M. M can be singular even though norm(M) ~=0. Of course, norm(M) == 0 means M is singular, but that condition seems overly restrictive (I think norm(M) == 0 is only true for the zero matrix?) .

M = [1 0;0 0];
norm(M)

ans = 1

det(M)

ans = 0

Am I missing something?

[Bruno Luong]

And of course det(M(w)) can converge to 0, "independent" of smalest/largest sinigular value

M1(w) = w*eye(n)

and

M2(w) = diag(1/w, w)

The question as it stands has very little interest mathematically, and if M is numercial approximated, one can never estimate accuracy smalleste eigen value below eps(1)*norm(M); which prevent ANY numerical algorithm to solve the problem OP asks.

In other word asking question on det is mostly a dead end numericaly in gerenal.

[Matt J]

I've edited my post to recommend rcond(M) rather than norm(M). As Bruno's examples suggest, rcond would also be problematic without an upper bound on the largest eigenvalue. However, from the OP's code, it appears that the M(i,j) in this case are all globally bounded functions of w. They are just sums of complex exponential terms whos phases are functions of w. That would be enough to globally bound the largest eigenvalue.

[Bruno Luong]

On a side note, numerically it is more reasonable to work with log(det(A)). It has some physical interpretation in some circumstance and equivalent to det(log(A)) such as A is definite positive matrix.

If the eigen values are purely complex phase, log(det) transforms to the sim of the phases.

[ma Jack]

Thanks for the discussion, I have learned a lot and I would like to add a fact here which might help to solve my problem better:

My goal is to find a complex w such that det(ModeXY)->0 (since this ensures that the eigenvector of an eigenequation like H\psi=lam\psi is not a 0-vector), but note that the w found is in the range where its real part is around 1 (roughly 0.6 to 1.5) and its imaginary part is negative (roughly -0.1 to 0), since my problem is relevant to the physical world, the correct w must be within this range, based on experimental observations.

[Bruno Luong]

What is the unit and physical interpretation of the determinant of M for you?

If you can't answer my question accurately I would bet you stil not totally clear about the claim that you know how to translate the physics into maths.

Ifthe unit is (eV)^8 you might also have trouble and think about the physics again.

The only possibility of valid answer is det(M) has unitless unit. Is it the case?

[ma Jack]

Uh .... I'm sorry, I didn't think too hard about this, but these formulas should be fine, I've seen people looking for complex roots on these formulas, but I'm not sure what method they used to find them.

[Matt J]

My goal is to find a complex w such that det(ModeXY)->0 (since this ensures that the eigenvector of an eigenequation like H\psi=lam\psi is not a 0-vector)

Eigenvectors are never 0-vectors, regardless of the determinant of H.

If you meant that you want to make sure H has a trivial null-space, det(H)=0 would be the opposite of what you want.

[ma Jack]

Sorry I did't make it clear, not det(H), but det(H-eye*lam), in addition I use norm found that the effect is not det good (I scaled the problem so that det can work).

[Matt J]

Sorry I did't make it clear, not det(H), but det(H-eye*lam),

Still not clear. The condition det(H-eye*lam)=0 ensures that lam is an eigenvalue. Is that what you want?

in addition I use norm found that the effect is not det good

Of course not. Above, we said the idea of using norm(M) should be abandoned. What about rcond?

[Bruno Luong]

det(H-eye*lam) (= 0)

So your problem is find the eigen values/vectors or H. Sorry but why not simply tell us that at the first place instead of asking a question with innacurate description that is not mathematically equivalent?

You know eig function in symbolic? or a numerical version of it?

[ma Jack]

Sorry please give me some time to redescribe the problem again, also the above code does not reflect my problem very accurately (due to simplification), I am now testing the difference between the various methods, it will take a while.

[Bruno Luong]

"due to simplification"

Please never do that when asking question. You would think it helps? Generally no, you should ask accurately question, since a simplified problem can have completily solution.

[ma Jack]

Hi everyone, interested in following my latest issue?I encountered a very, very strange problem when using the fminsearch+rcond function(Same code, same computer, same matlab version, different running times and different results).

[Rosalinda]

hello ma jack..i encounter the same problem for 8by 8 matrix..if you could please tell me how you resolve your issue

[Torsten]

Since @ma Jack 's problem description was very poor until the end, maybe you could again describe what you consider as "the same problem".