# Struggling to Vectorize the following Code.

3 views (last 30 days)
Jake Simmonds on 9 Nov 2019
Edited: Fabio Freschi on 9 Nov 2019
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

Fabio Freschi on 9 Nov 2019
Edited: Fabio Freschi on 9 Nov 2019
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

Fabio Freschi on 9 Nov 2019
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 on 9 Nov 2019
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 on 9 Nov 2019
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

Walter Roberson on 9 Nov 2019
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 on 9 Nov 2019
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 on 9 Nov 2019
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 on 9 Nov 2019
The bottleneck of the code is the plot. In this sense, plotting all data in a single call can be seen as vectorization.