fast evaluation of x(k+1) = x(k)*e(k)

Is there any special trick how to evaluate fast this recurrent formula:

x(1) = x0
x(k+1) = x(k)*e(k), k= 1,2, ...N-1 

where N = length(x) = length(e)+1

e ... known vector

I mean faster then standard for-loop:

x(1) = x0;
for k = 1:N-1
    x(k+1) = x(k)*e(k);
end

[Matt J]

If N is large enough that speed would matter, you are likely running the risk of underflow or overflow of the successive products.

[Michal]

@Matt J Underflow or overflow is not an issue in my case, bacause 0 <= e(k) <= 1

Answer 1

[Dyuman Joshi]

Timings of for loop with preallocation and vectorization are comparable for a large number of elements -

x = 69420;
N = 5e6;
e = rand(1,N);

F1 = @() forloop(x, e, N);
F2 = @() forlooppre(x, e, N);
F3 = @() vectorization(x, e);

fprintf('Time taken by for loop without preallocation for a vector with %g elements = %f seconds', N, timeit(F1))

Time taken by for loop without preallocation for a vector with 5e+06 elements = 0.320160 seconds

fprintf('Time taken by for loop with preallocation approach for a vector with %g elements = %f seconds', N, timeit(F2))

Time taken by for loop with preallocation approach for a vector with 5e+06 elements = 0.034487 seconds

fprintf('Time taken by vectorized approach for a vector with %g elements = %f seconds', N, timeit(F3))

Time taken by vectorized approach for a vector with 5e+06 elements = 0.036738 seconds

y1 = F1();
y2 = F2();
y3 = F3();

%Comparing results using tolerance
all(abs(y1-y2)<1e-10 & abs(y2-y3)<1e-10)

ans = logical
   1

%% Function definitions

%For loop without preallocation
function x = forloop(x, e, N)
    for k=1:N
        x(k+1) = x(k)*e(k);
    end
end

%For loop without preallocation
function x = forlooppre(x, e, N)
    %preallocation
    x = [x zeros(1, N)];
    for k=1:N-1
        x(k+1) = x(k)*e(k);
    end
end
    
%Vectorized approach
function x = vectorization(x, e)
    x = x.*[1 cumprod(e)];
end

[Michal]

Of course, I used the preallocation in my for-loop I just forgot to add it to my question :). Looks like, that for-loop with preallocation is still a bit faster (even for larga data) than fully vectorized approach, which is a bit surprising, too. Probably due to the cumprod overhead.

[Matt J]

Probably due to the cumprod overhead.

No, probably due to the overhead of horzcat.

[Dyuman Joshi]

I brought the concatenation out of the function and ran the timings, it looks like horzcat is the culprit for overhead -

x = 69420;
N = 5e6;
e = rand(1,N);
E = [1 e];
F1 = @() forloop(x, e, N);
F2 = @() forlooppre(x, e, N);
F3 = @() vectorization(x, E);

fprintf('Time taken by for loop without preallocation for a vector with %g elements = %f seconds', N, timeit(F1))

Time taken by for loop without preallocation for a vector with 5e+06 elements = 0.321753 seconds

fprintf('Time taken by for loop with preallocation for a vector with %g elements = %f seconds', N, timeit(F2))

Time taken by for loop with preallocation for a vector with 5e+06 elements = 0.035463 seconds

fprintf('Time taken by vectorized approach for a vector with %g elements = %f seconds', N, timeit(F3))

Time taken by vectorized approach for a vector with 5e+06 elements = 0.020578 seconds

y1 = F1();
y2 = F2();
y3 = F3();

%Comparing results using tolerance
all(abs(y1-y2)<1e-10 & abs(y2-y3)<1e-10)

ans = logical
   1

%% Function definitions

%For loop without preallocation
function x = forloop(x, e, N)
    for k=1:N
        x(k+1) = x(k)*e(k);
    end
end

%For loop without preallocation
function x = forlooppre(x, e, N)
    %preallocation
    x = [x zeros(1, N)];
    for k=1:N-1
        x(k+1) = x(k)*e(k);
    end
end
    
%Vectorized approach
function y = vectorization(x, e)
    y = x.*cumprod(e);
end

[Michal]

@Dyuman Joshi Thanks, your vectorized approach is clear winner now :)

[Matt J]

@Michal please accept-click the answer to indicate your question is resolved.

[Michal]

There is one remaining problem:

try to change:

x = 6942000

In this case are results different for tolerance 1e-10. So there is different numerical error propagation behavior between for-loop and vectorized approach. Any idea how to eliminate this problem?

[Dyuman Joshi]

Seems to be working correctly here -

x = 6942000;
N = 5e6;
e = rand(1,N);
E = [1 e];
F1 = @() forloop(x, e, N);
F2 = @() forlooppre(x, e, N);
F3 = @() vectorization(x, E);

fprintf('Time taken by for loop without preallocation for a vector with %g elements = %f seconds', N, timeit(F1))

Time taken by for loop without preallocation for a vector with 5e+06 elements = 0.316304 seconds

fprintf('Time taken by for loop with preallocation for a vector with %g elements = %f seconds', N, timeit(F2))

Time taken by for loop with preallocation for a vector with 5e+06 elements = 0.034462 seconds

fprintf('Time taken by vectorized approach for a vector with %g elements = %f seconds', N, timeit(F3))

Time taken by vectorized approach for a vector with 5e+06 elements = 0.010307 seconds

y1 = F1();
y2 = F2();
y3 = F3();
%Comparing results using tolerance
all(abs(y1-y2)<1e-10 & abs(y2-y3)<1e-10)

ans = logical
   1

%% Function definitions

%For loop without preallocation
function x = forloop(x, e, N)
    for k=1:N
        x(k+1) = x(k)*e(k);
    end
end

%For loop without preallocation
function x = forlooppre(x, e, N)
    %preallocation
    x = [x zeros(1, N)];
    for k=1:N-1
        x(k+1) = x(k)*e(k);
    end
end
    
%Vectorized approach
function x = vectorization(x, e)
    x = x.*cumprod(e);
end

[Michal]

What OS and Matlab release you are using?

[Dyuman Joshi]

The code above was ran on

ver

-----------------------------------------------------------------------------------------------------
MATLAB Version: 23.2.0.2463057 (R2023b) Update 5
MATLAB License Number: 0
Operating System: Linux 5.4.233-0504233-generic #202302250651 SMP Sat Feb 25 12:26:27 UTC 2023 x86_64
Java Version: Java 1.8.0_292-b10 with AdoptOpenJDK OpenJDK 64-Bit Server VM mixed mode
-----------------------------------------------------------------------------------------------------
MATLAB                                                Version 23.2        (R2023b)
Simulink                                              Version 23.2        (R2023b)
5G Toolbox                                            Version 23.2        (R2023b)
Aerospace Blockset                                    Version 23.2        (R2023b)
Aerospace Toolbox                                     Version 23.2        (R2023b)
Antenna Toolbox                                       Version 23.2        (R2023b)
Audio Toolbox                                         Version 23.2        (R2023b)
Automated Driving Toolbox                             Version 23.2        (R2023b)
Bioinformatics Toolbox                                Version 23.2        (R2023b)
Bluetooth Toolbox                                     Version 23.2        (R2023b)
Communications Toolbox                                Version 23.2        (R2023b)
Computer Vision Toolbox                               Version 23.2        (R2023b)
Control System Toolbox                                Version 23.2        (R2023b)
Curve Fitting Toolbox                                 Version 23.2        (R2023b)
DO Qualification Kit                                  Version 23.2        (R2023b)
DSP System Toolbox                                    Version 23.2        (R2023b)
Database Toolbox                                      Version 23.2        (R2023b)
Datafeed Toolbox                                      Version 23.2        (R2023b)
Deep Learning Toolbox                                 Version 23.2        (R2023b)
Econometrics Toolbox                                  Version 23.2        (R2023b)
Embedded Coder                                        Version 23.2        (R2023b)
Filter Design HDL Coder                               Version 23.2        (R2023b)
Financial Instruments Toolbox                         Version 23.2        (R2023b)
Financial Toolbox                                     Version 23.2        (R2023b)
Fixed-Point Designer                                  Version 23.2        (R2023b)
Fuzzy Logic Toolbox                                   Version 23.2        (R2023b)
Global Optimization Toolbox                           Version 23.2        (R2023b)
HDL Coder                                             Version 23.2        (R2023b)
HDL Verifier                                          Version 23.2        (R2023b)
IEC Certification Kit                                 Version 23.2        (R2023b)
Image Acquisition Toolbox                             Version 23.2        (R2023b)
Image Processing Toolbox                              Version 23.2        (R2023b)
Industrial Communication Toolbox                      Version 23.2        (R2023b)
Instrument Control Toolbox                            Version 23.2        (R2023b)
LTE Toolbox                                           Version 23.2        (R2023b)
MATLAB Compiler                                       Version 23.2        (R2023b)
MATLAB Compiler SDK                                   Version 23.2        (R2023b)
MATLAB Report Generator                               Version 23.2        (R2023b)
Mapping Toolbox                                       Version 23.2        (R2023b)
Mixed-Signal Blockset                                 Version 23.2        (R2023b)
Model Predictive Control Toolbox                      Version 23.2        (R2023b)
Navigation Toolbox                                    Version 23.2        (R2023b)
Optimization Toolbox                                  Version 23.2        (R2023b)
Parallel Computing Toolbox                            Version 23.2        (R2023b)
Partial Differential Equation Toolbox                 Version 23.2        (R2023b)
Phased Array System Toolbox                           Version 23.2        (R2023b)
Powertrain Blockset                                   Version 23.2        (R2023b)
Predictive Maintenance Toolbox                        Version 23.2        (R2023b)
RF Blockset                                           Version 23.2        (R2023b)
RF Toolbox                                            Version 23.2        (R2023b)
Requirements Toolbox                                  Version 23.2        (R2023b)
Risk Management Toolbox                               Version 23.2        (R2023b)
Robotics System Toolbox                               Version 23.2        (R2023b)
Robust Control Toolbox                                Version 23.2        (R2023b)
Sensor Fusion and Tracking Toolbox                    Version 23.2        (R2023b)
SerDes Toolbox                                        Version 23.2        (R2023b)
Signal Processing Toolbox                             Version 23.2        (R2023b)
SimBiology                                            Version 23.2        (R2023b)
SimEvents                                             Version 23.2        (R2023b)
Simscape                                              Version 23.2        (R2023b)
Simscape Driveline                                    Version 23.2        (R2023b)
Simscape Electrical                                   Version 23.2        (R2023b)
Simscape Fluids                                       Version 23.2        (R2023b)
Simscape Multibody                                    Version 23.2        (R2023b)
Simulink 3D Animation                                 Version 23.2        (R2023b)
Simulink Check                                        Version 23.2        (R2023b)
Simulink Code Inspector                               Version 23.2        (R2023b)
Simulink Coder                                        Version 23.2        (R2023b)
Simulink Control Design                               Version 23.2        (R2023b)
Simulink Coverage                                     Version 23.2        (R2023b)
Simulink Design Optimization                          Version 23.2        (R2023b)
Simulink Design Verifier                              Version 23.2        (R2023b)
Simulink Desktop Real-Time                            Version 23.2        (R2023b)
Simulink PLC Coder                                    Version 23.2        (R2023b)
Simulink Real-Time                                    Version 23.2        (R2023b)
Simulink Report Generator                             Version 23.2        (R2023b)
Simulink Test                                         Version 23.2        (R2023b)
Stateflow                                             Version 23.2        (R2023b)
Statistics and Machine Learning Toolbox               Version 23.2        (R2023b)
Symbolic Math Toolbox                                 Version 23.2        (R2023b)
System Identification Toolbox                         Version 23.2        (R2023b)
Text Analytics Toolbox                                Version 23.2        (R2023b)
Vehicle Dynamics Blockset                             Version 23.2        (R2023b)
Vehicle Network Toolbox                               Version 23.2        (R2023b)
Vision HDL Toolbox                                    Version 23.2        (R2023b)
WLAN Toolbox                                          Version 23.2        (R2023b)
Wavelet Toolbox                                       Version 23.2        (R2023b)
Wireless HDL Toolbox                                  Version 23.2        (R2023b)

Linux and MATLAB Version R2023b Update 5

[Michal]

Very strange!? The same code on R2023b_upd5 and my both machines LinuxMInt 21.2 +Windows 11 Pro machine shows:

for x = 6942000;

Time taken by for loop without preallocation for a vector with 5e+06 elements = 0.344369 seconds 
Time taken by for loop with preallocation for a vector with 5e+06 elements = 0.025052 seconds 
Time taken by vectorized approach for a vector with 5e+06 elements = 0.014689 seconds 

ans =

logical

0 

and

>> max( abs(y1-y2))
ans =
     0
     
>> max( abs(y1-y3))
ans =
4.6566e-10

for x = 69420;

are tolerance checks OK

Just try more independent runs with different random vector e!!!