Struggling to Vectorize the following Code.

This code is running really slow would like to vectorize it, however i don't really know how to do so in this case. The code is below, it does work but the N value really should be in the 1000's I've managed to vectorize code before for a much simplier system but this one i have no idea :( any help or suggestions are appreciated. Thanks

hold on

% Number of iterations
N = 100;
x = 1; y = 1;z = 1;
for w21 = linspace(-12,-3,N)
for i = 1:N-1
y = y_iterate(x,z,w21);
z = z_iterate(y);  
x = x_iterate(y);   
if i >= (N - 200)   
p = plot(w21,x,'.k','MarkerSize',3);  
        
end
    
end

end

function val = x_iterate(y)  
val = -3 + 8.*(1 ./ (1 + exp(-y))); 
end


function val = z_iterate(y) 
val = -7 + 8.*(1 ./ (1 + exp(-y)));
end

function val = y_iterate(x,z,w21)
val = 4 + w21.*(1 ./ (1 + exp(-x))) + 6.*(1 ./ (1 + exp(-z)));
end

Answer 1

[Fabio Freschi]

Almost all execution time is due to the plot function as highlighted by the profiler:

You can call the plot only once at the end of the two loops. One way is to load the data to plot in two vectors and call plot only once

hold on
% Number of iterations
N = 100;
x = 1; y = 1;z = 1;

% initialization
aa = [];
bb = [];
for w21 = linspace(-12,-3,N)
    for i = 1:N-1
        y = y_iterate(x,z,w21);
        z = z_iterate(y);
        x = x_iterate(y);
        if i >= (N - 200)
            % update
            aa = [aa; w21];
            bb = [bb; x];            
        end      
    end
end
p = plot(aa,bb,'.k','MarkerSize',3);

function val = x_iterate(y)
val = -3 + 8.*(1 ./ (1 + exp(-y)));
end
function val = z_iterate(y)
val = -7 + 8.*(1 ./ (1 + exp(-y)));
end
function val = y_iterate(x,z,w21)
val = 4 + w21.*(1 ./ (1 + exp(-x))) + 6.*(1 ./ (1 + exp(-z)));
end

[Jake Simmonds]

Hi thanks for responding it still seems to run really slow even with plotting just once! Not sure how to speed it up i was recommended using the following:

[X,Y,Z] = deal( zeros(N,N-1) );

w21 = linspace(-12,-3,N)';

for i = 1:N-1
    
y = y_iterate(x,z,w21);
    
z = z_iterate(y);
    
x = x_iterate(y);
    
X(:,i) = x;
    
Y(:,i) = y;
    
Z(:,i) = z;

end

plot(w21,X)

[Jake Simmonds]

I think the above plots incorrectly though :(

[Fabio Freschi]

Sorry, the plot function it is mistyped, the correct version is

p = plot(aa,bb,'.k','MarkerSize',3);

Regarding the timing, it is faster:

Old: Elapsed time is 7.446350 seconds.

New: Elapsed time is 0.169045 seconds.

[Jake Simmonds]

Using p = plot(aa,bb,'.k','MarkerSize',3); it still runs quite slow ! the above code in my comment runs really quick but it seems to plot to many points for some reason.

[Walter Roberson]

[X,Y,Z] = deal( zeros(N,N-1) );

That should probably be

[X,Y,Z] = deal( zeros(1,N-1) );

unless your initial x or z are vectors. If they are vectors then

[X,Y,Z] = deal( zeros(length(x),N-1) );

[Fabio Freschi]

I see very similar plots from your code and my code

Old

New

and the timings are very different:

Old: Elapsed time is 7.446350 seconds.

New: Elapsed time is 0.169045 seconds.

[Fabio Freschi]

To speed up further, can you guess the dimension of the vectors to plot as a function of the parameter N? we can improve the performances also for large N.

[Fabio Freschi]

Try this: for N = 1000 it takes about 4 secs

hold on
% Number of iterations
N = 1000;
x = 1; y = 1;z = 1;
% initialization
aa = zeros(200*N,1);
bb = zeros(200*N,1);
k = 1;

for w21 = linspace(-12,-3,N)
    for i = 1:N-1
        y = y_iterate(x,z,w21);
        z = z_iterate(y);
        x = x_iterate(y);
        if length(x) ~= 1
            p = 1;
        end
        if i >= (N - 200)
            % update
            aa(k) = w21;
            bb(k) = x;
            k = k+1;
        end      
    end
end
p = plot(aa(1:k-1),bb(1:k-1),'.k','MarkerSize',3);

function val = x_iterate(y)
val = -3 + 8.*(1 ./ (1 + exp(-y)));
end
function val = z_iterate(y)
val = -7 + 8.*(1 ./ (1 + exp(-y)));
end
function val = y_iterate(x,z,w21)
val = 4 + w21.*(1 ./ (1 + exp(-x))) + 6.*(1 ./ (1 + exp(-z)));
end

Answer 2

[Walter Roberson]

No, you have a chaotic system. Each iteration depends on the bitwise details of the previous iteration.

It is possible to calculate the algebraic iteration to any given number of steps in terms of the original inputs, getting out more complicated formulas each step, but the rounding that takes place at each individual step turns out to be important to the long-term outcome, and the algebraic formulas for multiple steps will not reproduce that round-off the same way.

Example: after two iterations, the new y is

w21/(exp(3 - 8/(exp(- w21/(exp(-x) + 1) - 6/(exp(-z) + 1) - 4) + 1)) + 1) + 6/(exp(3 - 8/(exp(- w21/(exp(-x) + 1) - 6/(exp(-z) + 1) - 4) + 1)) + 1) + 4

where those are the initial x, z values in the formula.

Your system immediately loses its original y coordinate, by the way: the original y coordinate plays no role in the system.

[Jake Simmonds]

Not sure what your point is here i'm afraid? I'm plotting a bifurication diagram for a dynamical system i understand the mathematics of what i'm doing , i just don't understand how to code it so that it plots the diagram efficently and not slowly.

[Walter Roberson]

Vectorization would be the calculation of a number of steps at the same time, possibly in parallel behind the scenes, the way that sin(x) for a large vector of x can be split up between cores and calculated by hard-coded routines and the results put back together.

You cannot do that kind of vectorization because every step depends upon the bit-for-bit details of every other step.

The closest you might be able to get to vectorization might be working on the GPU. See for example https://www.mathworks.com/company/newsletters/articles/gpu-enables-obsession-with-fractals.html

[Fabio Freschi]

The bottleneck of the code is the plot. In this sense, plotting all data in a single call can be seen as vectorization.