**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

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

2 views (last 30 days)

Show older comments

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.

##### 6 Comments

ma Jack
on 6 Jul 2022

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

ma Jack
on 6 Jul 2022

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
on 6 Jul 2022

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.

### Answers (1)

Matt J
on 6 Jul 2022

Edited: Matt J
on 7 Jul 2022

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)

##### 19 Comments

ma Jack
on 6 Jul 2022

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
on 6 Jul 2022

Edited: Matt J
on 6 Jul 2022

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
on 6 Jul 2022

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
on 6 Jul 2022

Paul
on 7 Jul 2022

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
on 7 Jul 2022

Edited: Bruno Luong
on 7 Jul 2022

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
on 7 Jul 2022

Edited: Matt J
on 7 Jul 2022

Bruno Luong
on 7 Jul 2022

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
on 7 Jul 2022

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
on 7 Jul 2022

Edited: Bruno Luong
on 7 Jul 2022

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
on 7 Jul 2022

Matt J
on 7 Jul 2022

Edited: Matt J
on 7 Jul 2022

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
on 8 Jul 2022

Matt J
on 8 Jul 2022

Edited: Matt J
on 8 Jul 2022

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
on 8 Jul 2022

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?

ma Jack
on 8 Jul 2022

Bruno Luong
on 8 Jul 2022

"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
on 12 Jul 2022

### See Also

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

### Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)