Why do we get different results here?

I have two functions namely "myfunAskMathworks.m" and "ByAskic.m". The latter is a modified version of the 1st one. In the 1st one, the antenna array is constructed within these lines:

vecH = @(MAT) MAT(:).'; 
steerVecT = @(ang) exp(1j*2*pi*d*(0:M-1).'*sin(vecH(ang))); 
steerVecR = @(ang) exp(1j*2*pi*d*(0:N-1).'*sin(vecH(ang)));

while in the 2nd code the antenna array is constructed separately by another m file called "ULA.m" which is called inside the 2nd code i.e., within ByAskic.m.

In addition to these two m files, I have three other supporting m-files namely "fpa1.m", "main.m" and "ULA.m".

PROBLEM: The problem is that when I run the m-file "main.m" and check my results in command window, the results obtained with the 1st function "myfunAskMathworks.m" are correct but the results obtained with the 2nd function "ByAskic.m" are wrong.

How to check Results: Inside "main.m", u is my desired vector and I want to get this vector exactly or nearly like this also in the result after running the main.m. The result will be returned in a variable called 'two1". When I run the 'main.m', with the function "myfunAskMathworks.m", my result is correct but when I run the same 'main.m' with the function "ByAskic.m", my result is wrong.

How to run main.m: Inside main.m, you will see these two lines on line 15 and line 16 respectively:

   [best,fmin,time]=fpa1(10,0.8,2000,dim,lb,ub,@(b)ByAskic(b,u));
  %[best,fmin,time]=fpa1(10,0.8,2000,dim,lb,ub,@(b)myfunAskMathworks(b,u));

When the line 15 is uncommented like above, it means the 2nd function "ByAskic.m" is in use but when line 15 is commenred and line 16 is uncommented, then the 1st function "myfunAskMathworks.m" is in use.

Why it is so?

All the required files are in the attachment.

Answer 1

[Askic V]

Can you try with the following ByAskic.m function:

function e = ByAskic(b,u)
% u = [10 20 30 40];b = u;
N = 10;%6; % Tx antennas 
M = 10;% Rx antennas

rT=ULA(N);% Call ULA Array at Tx side
rR=ULA(M);% Call ULA Array at Rx side

% rT is size Nx3, so to make multiplication work, 
% [cos(alpha).*cos(beta); sin(alpha).*cos(beta); sin(beta)] must be 3x1
AT = @(alpha, beta) exp(-1j.*rT*pi*[cos(alpha).*cos(beta); sin(alpha).*cos(beta); sin(beta)]);

% rR is size Mx3, so to make multiplication work, 
% [cos(alpha).*cos(beta); sin(alpha).*cos(beta); sin(beta)] must be 3x1
AR = @(alpha, beta) exp(-1j.*rR*pi*[cos(alpha).*cos(beta); sin(alpha).*cos(beta); sin(beta)]);


%%%%%%%%%%%%%%%%%%%%
% Swapping vector b
%%%%%%%%%%%%%%%%%%%%
[~, ix] = sort(u); 
[~, ix1(ix)] = sort(b);  
b = b(ix1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculation of yo and ye
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yo = yMatTR(deg2rad(u), AT, AR); 
ye = yMatTR(deg2rad(b), AT, AR); 

%%%%%%%%%%%%%%%%%%
% MSE
%%%%%%%%%%%%%%%%%%
e = norm(yo-ye).^2/(M*N);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function yMatTR
%%%%%%%%%%%%%%%%%%%%%%%%5%%%
function y = yMatTR(targetAngle, steerVecT, steerVecR)
    az = deg2rad(targetAngle);
    el = deg2rad(targetAngle);
    steerA = steerVecT(az,el);
    steerB =  steerVecR(az,el);
    y=sum( steerA.*permute(steerB,[3,2,1]) ,2);
   y=y(:);   
end

[Sadiq Akbar]

Thanks a lot dear Askic V for your kind response. I tried your above code but it also gives wrong results.

[Askic V]

I really don't know what is the expected result or what is right or wrong result.

For example the following main function:

clear;clc;
u=[-55 55 65 -65];% Desired vector
dim=length(u);
lb=-90*ones(1,dim);
ub=90*ones(1,dim);

rng(0) % to get repeatable results
Runs=10;
one=zeros(Runs,1);
time1=zeros(Runs,1);
temp=zeros(Runs,dim);
two=zeros(Runs,dim);
nn=0;
for n=1:Runs %------------(2)
    nn=nn+1;
  [best,fmin,time]=fpa1(10,0.8,2000,dim,lb,ub,@(b)ByAskic(b,u));
  %[best,fmin,time]=fpa1(10,0.8,2000,dim,lb,ub,@(b)myfunAskMathworks(b,u));
  
  one(nn)=fmin;
  temp(nn,:)=best;
  time1(nn)=time; 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Swapping
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  [~, ix] = sort(u);
  [~, ix1(ix)] = sort(temp(nn,:));  
  two(nn,:) = temp(nn,ix1); 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Save the Workspace data 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%save myData  
%%%%%%%%%%%%%%%%%%%%%%%%
% Checking Results in command window
%%%%%%%%%%%%%%%%%%%%%%%%
[one1 ind]=sort(one,'descend');
[fitness,ind1]=min(one1)
two1=two(ind1,:)
u

with the byAsikc function:

function e=ByAskic(b,u)
% u = [10 20 30 40];b = u;
N = 10;%6; % Tx antennas 
M = 10;% Rx antennas

rT=ULA(N);% Call ULA Array at Tx side
rR=ULA(M);% Call ULA Array at Rx side

% rT is size Nx3, so to make multiplication work, 
% [cos(alpha).*cos(beta); sin(alpha).*cos(beta); sin(beta)] must be 3x1
AT = @(alpha, beta) exp(-1j.*rT*pi*[cos(alpha).*cos(beta); sin(alpha).*cos(beta); sin(beta)]);

% rR is size Mx3, so to make multiplication work, 
% [cos(alpha).*cos(beta); sin(alpha).*cos(beta); sin(beta)] must be 3x1
AR = @(alpha, beta) exp(-1j.*rR*pi*[cos(alpha).*cos(beta); sin(alpha).*cos(beta); sin(beta)]);
%%%%%%%%%%%%%%%%%%%%
% Swapping vector b
%%%%%%%%%%%%%%%%%%%%
[~, ix] = sort(u); 
[~, ix1(ix)] = sort(b);  
b = b(ix1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculation of yo and ye
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yo = yMatTR(u, AT, AR); 
ye = yMatTR(b, AT, AR); 
%%%%%%%%%%%%%%%%%%
% MSE
%%%%%%%%%%%%%%%%%%
e=norm(yo-ye).^2/(M*N);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function yMatTR
%%%%%%%%%%%%%%%%%%%%%%%%5%%%
function y = yMatTR(targetAngle, steerVecT, steerVecR)
    az = deg2rad(targetAngle);
    el = deg2rad(targetAngle);
    steerA = steerVecT(az,el);
    steerB =  steerVecR(az,el);
    y=sum( steerA.*permute(steerB,[3,2,1]) ,2);
   y=y(:);   
end

In my case (note added rng(0) ) will produce the following output:

two1 =

 -54.720454966995689  55.336380456191215  65.757549351373029 -64.346201757620520

u =
   -55    55    65   -65

But calling the other function:

    [best,fmin,time]=fpa1(10,0.8,2000,dim,lb,ub,@(b)myfunAskMathworks(b,u));

produce the following result:

    two1 =
 -25.016307607428871  35.004875903879281  65.001676522539057 -35.023254511076644
u =
   -55    55    65   -65

So in my case, the bigger deviations are with myfunAsk... function

[Sadiq Akbar]

@Askic V

Thanks a lot dear Askic V for your kind response and sorry for my late response as I was out of station. I copied your above code and ran it several times but each time it gives me the following result:

two1 = -64.9891 -54.9968 55.0031 -65.0108

u = -55 55 65 -65

If you see my desired result should be like u or nearly like u but I get the desired result in two1 which is not like u or nearly like u.

Though I run your code several times but each time I get the above result which is not like u or nearly like u.

[Askic V]

You're using essentially optimization algorithm that has random part built in. With different seeds, the results are indeed different. What are the result if you call myfunAskMathworks(b,u)?

Did you put rng(0) just before Runs = 10; line?

Anyway, I have noticed that both functions when called in fpa1 can return very different results.

I'm not familiar with fpa1 function, and I think I canno help you further.

I have run this once again and obtained follwoing results:

% for ByAskic(b,u)
two1 =
 -54.720454966995689  55.336380456191215  65.757549351373029 -64.346201757620520

u =
   -55    55    65   -65
   
% [best,fmin,time]=fpa1(10,0.8,2000,dim,lb,ub,@(b)ByAskic(b,u));
two1 =
 -25.016307607428871  35.004875903879281  65.001676522539057 -35.023254511076644

u =
   -55    55    65   -65

[Sadiq Akbar]

Thanks a lot dear Askic V for your kind response.

What are the result if you call myfunAskMathworks(b,u)? I ran it several times but it gives the correct result i.e., it gives the desired result i.e.,

two1 = -54.8890 54.6959 64.7207 -64.9280

u = -55 55 65 -65

But when I run it with "ByAskic.m", it gives me wrong result though the same optimization algorithm is used that has random part built in. I don't know why?

Also if wrong answers is caused by fpa1 i.e., the optimization algorithm that has random part built in then the wrong answers should be obtained in both the cases but it's not so here.

Answer 2

[Askic V]

@Sadiq Akbar I'm puzzled too. In general, the results are not correct event sometimes are.

For example, the result is wrong for input vector u = [10,20,30,40], but the result will be correct if lb is in that case [0 0 0 0].

I have reduced runs to 1 and increased psize from 10 to 25 and get results as I described above (didn't attached updated main.m). That is why I put all functions to the one script.

This is the oputput of the following script:

two1 = 1×4

-54.4730 55.5201 65.8900 -64.1584

clear;clc;
u=[-55, 55, 65, -65];% Desired vector
dim=length(u);
lb=-90*ones(1,dim);
ub=90*ones(1,dim);

rng(0) % to get repeatable results
Runs=1;
one=zeros(Runs,1);
time1=zeros(Runs,1);
temp=zeros(Runs,dim);
two=zeros(Runs,dim);
nn=0;
for n=1:Runs %------------(2)
    nn=nn+1;
  [best,fmin,time]=fpa1(25,0.8,2000,dim,lb,ub,@(b)ByAskic(b,u));
  
  one(nn)=fmin;
  temp(nn,:)=best;
  time1(nn)=time; 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Swapping
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  [~, ix] = sort(u);
  [~, ix1(ix)] = sort(temp(nn,:));  
  two(nn,:) = temp(nn,ix1); 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Save the Workspace data 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%save myData  
%%%%%%%%%%%%%%%%%%%%%%%%
% Checking Results in command window
%%%%%%%%%%%%%%%%%%%%%%%%
[one1 ind]=sort(one,'descend');
[fitness,ind1]=min(one1)

fitness = 5.7560e-04

ind1 = 1

two1=two(ind1,:)

two1 = 1×4

-54.4730 55.5201 65.8900 -64.1584

u

u = 1×4

-55 55 65 -65

function [best,fmin,time]=fpa1(psize,psw,iter,dim,lb,ub,obj)
    tic;
    fun=obj;
    n=psize;           % Population size, typically 10 to 25
    p=psw;             % probabibility switch
    N_iter=iter;       % Total number of iterations
    d=dim;
    Lb=ones(1,d).*lb;  %Lower bounds.
    Ub=ones(1,d).*ub;  %Upper bounds 

    Sol=zeros(psize,dim); % suggested by Dyuman Joshi for speed up
    Fitness=zeros(psize,1);% suggested by Dyuman Joshi for speed up

    for i=1:n
        Sol(i,:)=Lb+(Ub-Lb).*rand(1,d);
        Fitness(i)=fun(Sol(i,:));%--------------------(1)  
    end

    [fmin,I]=min(Fitness);
    best=Sol(I,:);
    S=Sol;
    for t=1:N_iter
        for i=1:n
            if rand>p
               L=Levy(d);
               dS=L.*(Sol(i,:)-best);
               S(i,:)=Sol(i,:)+dS;
               S(i,:) = min(max(S(i,:), Lb), Ub);%S(i,:)=simplebounds(S(i,:),Lb,Ub);%-By Walter Roberson           
            else
                epsilon=rand;
                JK=randperm(n, 2);%JK=randperm(n); %----------------By Walter Roberson           
                S(i,:)=S(i,:)+epsilon*(Sol(JK(1),:)-Sol(JK(2),:));
                S(i,:) = min(max(S(i,:), Lb), Ub);%S(i,:)=simplebounds(S(i,:),Lb,Ub);%-By Walter Roberson
            end        
            Fnew=fun(S(i,:)); %--------------------(2)   
            if (Fnew<=Fitness(i))
                Sol(i,:)=S(i,:);
                Fitness(i)=Fnew;
            end        
            if Fnew<=fmin
                best=S(i,:)   ;
                fmin=Fnew   ;
            end
        end
        time=toc;
     end

end

function L=Levy(d)
    beta=3/2;
    sigma=(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);
    u=randn(1,d)*sigma;
    v=randn(1,d);
    step=u./abs(v).^(1/beta);
    L=0.01*step;
end

function e=ByAskic(b,u)
    % u = [10 20 30 40];b = u;
    N = 10;%6; % Tx antennas 
    M = 10;% Rx antennas

    rT=ULA(N);% Call ULA Array at Tx side
    rR=ULA(M);% Call ULA Array at Rx side

    % rT is size Nx3, so to make multiplication work, 
    % [cos(alpha).*cos(beta); sin(alpha).*cos(beta); sin(beta)] must be 3x1
    AT = @(alpha, beta) exp(-1j.*rT*pi*[cos(alpha).*cos(beta); sin(alpha).*cos(beta); sin(beta)]);

    % rR is size Mx3, so to make multiplication work, 
    % [cos(alpha).*cos(beta); sin(alpha).*cos(beta); sin(beta)] must be 3x1
    AR = @(alpha, beta) exp(-1j.*rR*pi*[cos(alpha).*cos(beta); sin(alpha).*cos(beta); sin(beta)]);
    %%%%%%%%%%%%%%%%%%%%
    % Swapping vector b
    %%%%%%%%%%%%%%%%%%%%
    [~, ix] = sort(u); 
    [~, ix1(ix)] = sort(b);  
    b = b(ix1);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Calculation of yo and ye
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    yo = yMatTR(u, AT, AR); 
    ye = yMatTR(b, AT, AR); 
    %%%%%%%%%%%%%%%%%%
    % MSE
    %%%%%%%%%%%%%%%%%%
    e=norm(yo-ye).^2/(M*N);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function yMatTR
%%%%%%%%%%%%%%%%%%%%%%%%5%%%
function y = yMatTR(targetAngle, steerVecT, steerVecR)
    az = deg2rad(targetAngle);
    el = deg2rad(targetAngle);
    steerA = steerVecT(az,el);
    steerB =  steerVecR(az,el);
    y=sum( steerA.*permute(steerB,[3,2,1]) ,2);
   y=y(:);   
end

function r=ULA(N)
    %N = 10; % Number of antennas
    r = [(-(N-1)/2:(N-1)/2).',zeros(N,2)]; % Assume uniform linear array
end