Questions about the Lsqnonlin function

2022/7/2

Hi all,

Are there any experts here who use the Lsqnonlin function?I'm trying to use the Lsqnonlin function to optimize a more complex function so that its result converges to 0 (i.e. I expect to find an input value that makes the function result most closely converge to 0),But it reports an error, it says error using eig function because the input matrix contains NAN or inf, I don't know how to fix it, can anyone give me an answer? thanks in advance!

By the way the function I want to optimize is a determinant, I found that as long as I reduce the size of the determinant to 4*4, Lsqnonlin does not report an error, but as soon as I make the size of the determinant bigger, it reports an error.

The following is the error information：

Wrong use of eig

The input matrix contains NaN or Inf.

error trust (line 33)

[V,D] = eig(H);

error trdog (line 132)

st = trust(rhs,MM,delta);

error snls

error lsqncommon (line 169)

snls(funfcn,xC,lb,ub,flags.verbosity,options,defaultopt,initVals.F,initVals.J,caller,

...

error lsqnonlin (line 258)

lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,allDefaultOpts,caller,...

2022/7/3

Hi all,

This is my code, it is a function file not a script file, where Zero_seek_Main is the main function, the other functions are sub-functions (in order to construct the matrix).

When I reduce the size of the GvXY matrix (4*4 and below), lsqnonlin does not report an error, which is a bit strange.

Looking forward to any reply.

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

2022/7/3 pm

Hi all,

Some people have asked me why I require the determinant of the matrix, and the explanation is as follows:

Obviously, it is a characteristic equation,But note that the eigenvector [a;b;c;d] is unknown, and omega is also unknown, in order to ensure that the eigenvector [a;b;c;d] is not [0;0;0;0], so I need to make det(Gv(w)-lam(w). *eye(4))=0, or converge to 0, because in practice I found it impossible to find the right omega to make det=0.And also I found that fsolve is not as effective as lsqnonlin, or not effective at all.

Dear experienced experts, how do you solve the above problems in your actual projects? It should be very common in practice, right?

[Torsten]

We can only repeat what the error message says:

During the optimization with lsqnonlin, you try to determine eigenvalues of a matrix which has NaN or Inf elements. At this stage, lsqnonlin stops and throws the error.

[ma Jack]

Hi Mr.Torsten

Thank you for your comments,I have also made efforts for this mistake.Before finding the determinant of a matrix I use the following operations on it:

Mat(isnan(Mat))=0;

Mat(isinf(Mat))=0;

This ensures that the matrix I construct has no NAN or INF,Then I find the determinant for it, substitute it into the Lsqnonlin function and it still reports an error (as described in the previous question).

I tried to make sense of the internal operations of the lsqnonlin function, but I gave up because the lsqnonlin function is too complex, but I can ensure that it is the internal calculations of the lsqnonlin function that generate the NAN or INF, and that the values I entered do not have a NAN or INF.

Expect a reply , and if necessary I'll post the relevant code here.

[ma Jack]

Is there a function like lsqnonlin in matlab? The problem I want to solve is a nonlinear problem. My goal is to find an independent variable that makes the output of the optimized function closest to 0. Note that the independent variable must be complex, which means that fzero cannot be used.

[Torsten]

The solver isn't the problem, I guess. But we cannot help without execuable code.

And by the way: the name of the solver is Lsqnonlin, not isqnonlin.

[Walter Roberson]

See https://www.mathworks.com/help/optim/ug/complex-numbers-in-optimization-toolbox-solvers.html

lsqnonlin() can handle complex values, under a number of restrictions.

However, if you are trying to find an independent value that makes the output of the optimized function closest to 0, then lnsqnonlin() is not appropropriate: lsqnonlin() is for fitting a model to a set of data. You should probably be using fsolve(), which is able to handle complex values under the restrictions described above.

There are other approaches, such as splitting real and complex parts into separate variables, and emitting norm() of the objective. Such techniques can permit bounds.

[Walter Roberson]

Finding a value that drives a determinent to 0... that sounds like an eigenvalue problem ?

[ma Jack]

Thanks for all the comments and answers,I was pleasantly surprised by the lively atmosphere here.Please give me some time to read your suggestions carefully and simplify my code (the original code is a bit complicated, about 500 lines), then I will post the code here.

[ma Jack]

Hi all,

My issue is updated and contains my code,looking forward to any suggestions.

[Torsten]

You are working with a matrix ModeXY with values in the order of 1e30 for both real and imaginary part of which you try to take the determinant. This must result in numerical nonsense.

Review your code for errors, scaling of the equations or a reformulation of your problem. This way at least, you won't succeed.

[ma Jack]

Hi Mr.Torsten,

Thanks for your approach, lsqnonlin is no longer reporting errors, but I still need to test extensively, does anyone else have any suggestions?

[Torsten]

Yes. Searching for w, [a,b,c,d] under the constaint a^2+b^2+c^2+d^2 = 1 should make the use of "det" superfluous.

[ma Jack]

Hi Mr.Torsten

Sir why your equation does not contain Gv(w) in it? It is also a function of w

[Torsten]

I don't understand what you mean.

Searching for a vector x=[a b c d] with norm 1 excludes the case of the trivial solution x = 0.

[ma Jack]

Sorry sir, I didn't describe the problem clearly, the vectors [a;b;c;d] are completely unknown (no expressions), only Gv(w) and lambda(w) are known (with specific formulas)

Answer 1

[John D'Errico]

This is NOT a problem with lsqnonlin, but in your understanding of determinants and why they are a bad thing to use.

You misunderstand the problem. You want to compute a determinant. Well, you THINK you want to do that. But determinants are a classically BAD thing to work with. They are numerically unstable for a variety of reasons. For example, do you agree that the determinant of

A = eye(1200);

has no infs or NaNs in it? Does it return something meaningful for the determinant?

det(A)

ans = 1

WOW!. It does. And that is exactly what I would have expected. Next, if we just multiply A by 2, does B have any infs or NaNs in it? (I hope not.)

B = 2*A;
det(B)

ans = Inf

Surely B had nothing in it that caused the problem. Could it have been lsqnonlin?

How about C?

C = A/2;
det(C)

ans = 0

Are any of those matrices singular, thus A, B, or C? Of course not. In fact, they are all extremely well behaved.

So, it must have been lsqnonlin that caused this problem? But, wait. I never used lsqnonlin here. I wonder what happened?

The problem in the above examples is that determinants have nasty scaling problems. But they are also bad things to use to use in general, even for smaller matrices. For example, consider this simple, smaller matrix.

D = rand(7,6)*1000;
D = D*rand(6,7);

I KNOW, ABSOLUTELY KNOW that this matrix is singular. It must have a zero determinant, right?

det(D)

ans = -4.9378e+03

And that was just a 7x7 matrix, composed of moderately small entries. Nothing too large. Nothing terribly small. But the determinat came out being pretty darn large, compared to the zero number you would be looking to find.

The point being that determinants are NOT a good way to decide if a matrix is singular. They are NOT a good thing to use for whatever your task may be.

Yes, I know, your teachers told you to use determinants. Their teachers told them the same thing. But nobody ever explained that in terms of numerical compuutations, the things are just useless, at least when double precision arithmetic is used to compute them. In fact, they are not just useless, but actively bad because they lead people to do silly things like this. That was the class your teachers probably forgot to give, even if they knew somewhere in the backs of their minds that sometimes there are problems. And of course, everything works great when you compute the determinant of a 3x3 integer matrix in that homework assignment for class. So what could possibly go wrong?

Don't believe me? Try this one, this time using integers. The nice thing is now I can use sym to compute the determinant too.

E = randi(10,[10,9])*randi(10,[9,10])

E = 10×10

262 200 271 311 399 261 412 424 239 373 179 161 222 248 275 220 269 260 227 289 216 182 216 238 240 177 282 227 150 271 267 222 299 307 392 298 431 408 300 380 280 188 299 283 378 255 412 391 259 377 137 139 170 236 226 131 232 272 148 290 195 164 242 275 261 187 294 288 225 336 319 246 349 361 458 334 460 436 310 442 198 178 239 289 346 225 330 346 260 341 232 235 258 291 345 256 325 294 255 344

Again, EXACT INTEGERS. But the way I created the matrix, I KNOW that it is singular. It has rank at most 9, for a 10x10 matrix since it is the product of 10x9 and 9x10 matrices.

Now compute the determinant in double precision arithmetic, and then using symbolic tools.

det(E)

ans = -26.3364

det(sym(E))

ans = 

0

Do you see the difference? You CANNOT effectively use a determinant for such a task, searching for a zero. They are poorly behaved, almost always.

Again, ABSOLUTELY NOTHING to do with lsqnonlin.

We can understand what is happening, if we look at an LU decomposition of these matrices.

[L,U] = lu(E);
format long g
diag(U)

ans = 10×1

319 56.0909090909091 27.6457273794221 42.1032597818726 -96.2713298870711 -30.455195763015 56.635421538508 -20.1280502645681 -29.7973765940712 -1.26970747187275e-14

Remember that L and U are lower and upper triangular matrices. As well, the diagonal elements of L are all exactly 1. So the determinant of E is just the product of the diagonal elements of U. Do you see that one of those numbers is not exactly zero, but only close to it. As close as we can get in double precision arithmetic.

prod(diag(U))

ans = 
         -26.3364250073676

So even though we know that E is singular, the product of those diagonal elements is pretty far away from zero. And so was the determinant.

Again, I don't know why you think you need to compute a determinant here. You don't want to do it. No matter what someone told you, or what you remember from some barely remembered class, you don't want to do it. Perhaps you can show what you are trying to do, then we could offer some advice on how to do it better. Lacking that, all I can say is DON'T DO WHAT YOU ARE DOING.

[ma Jack]

Thank you very much for your detailed answer sir, I can't believe that all these are true, it seems that my math knowledge is still too poor. Also I have updated my answer, it includes my code and why I asked for a value so that det(Matrix) = 0 (or tends to 0).