Vectorize MATLAB for loop

Question

I have the following lines of code

y = zeros(n, 1);

for i=1:n 
    b = L * [u(i:-1:max(1,i-M+1));zeros((-i+M)*(i-M<0),1)];
    y(i) = b' * gamma;
end

u is nx1, gamma is Mx1 and L is MxM

n takes very large values, so are there any ideas on how to vectorize the for loop?

please provide more details, what size is `u`? is `gamma` a scalar? — m.s., May 31 '15 at 11:30
then memory consumption can be an issue when trying to unroll the loop ... how long does the for loop take for you? — m.s., May 31 '15 at 12:21
I have written a version which is based on matrix multiplication without any loops, however due the memory consumption for your large `n` this does not help at all. you could try `parfor` or moving the calculation on to a GPU — m.s., May 31 '15 at 12:58

score 2 · Answer 1 · edited Jun 20 '20 at 09:12

Discussion and solution code

Initial approach

Matrix-multiplications based approach -

u_pad = [zeros(M-1,1) ; u];         %//   Pad u with zeros
idx = bsxfun(@plus,[M:-1:1]',0:n-1);%//'# Calculate sliding indices for u
u_pad_indexed = u_pad(idx);         %//   Index into padded u 
y_vectzed = gamma.'*L*u_pad_indexed;%//'# Matrix-multiplications for final o/p

Modified approach

Now, you have huge datasizes to work with. So, to optimize for such a case, the data could be broken into smaller pieces that are runnable and the operations could be done iteratively. Then, each iteration would calculate a part of the output array.

With this new strategy, the intial setting up could be done once and reused within each iteration. The modified approach would look something like this -

div_factor = [...]   %// Make sure it is a divisor of n
nrows = n/div_factor;
start_idx =  bsxfun(@plus,[M:-1:1]',0:nrows-1);  %//'
u_pad = [zeros(M-1,1) ; u];

y_vectorized = zeros(div_factor,n/div_factor);
for iter = 1:div_factor
    u_pad_indexed = u_pad((iter-1)*nrows + start_idx);
    y_vectorized(iter,:) = gamma.'*L*u_pad_indexed;  %//'
end
y_vectorized = reshape(y_vectorized.',[],1);

Benchmarking

%// Size parameters and input arrays
n = 4000000;
M = 1000;
u = rand(n,1);
gamma = rand(M,1);
L = rand(M,M);

%// Warm up tic/toc.
for k = 1:50000
    tic(); elapsed = toc();
end

disp('----------- With Original loopy code');
tic
y = zeros(n, 1);
for i=1:n
    b = L * [u(i:-1:max(1,i-M+1));zeros((-i+M)*(i-M<0),1)];
    y(i) = b' * gamma;  %//'
end
toc
clear b y

disp('----------- With Proposed solution code');
tic
..... Proposed Modified Code with div_factor = 200
toc

Runtimes

----------- With Original loopy code
Elapsed time is 498.563049 seconds.
----------- With Proposed solution code
Elapsed time is 44.273299 seconds.

user1543042 · Accepted Answer · 2015-05-31T18:51:12.770

Edit: I am redoing the solution because I found out that Matlab does not handle anonymous functions well. So I changed the call from an anonymous function to a normal function. Making this change:

Test 1

Comparison(40E3, 3E3)
Elapsed time is 21.731176 seconds.
Elapsed time is 251.327347 seconds.
|y2-y1| = 3.1519e-06

Test 2

Comparison(40E3, 1E3)
Elapsed time is 6.407259 seconds.
Elapsed time is 25.551116 seconds.
|y2-y1| = 2.8402e-07

Test 3

Comparison(20E3, 3E3)
Elapsed time is 10.484422 seconds.
Elapsed time is 125.033313 seconds.
|y2-y1| = 1.5462e-06

Test 4

Comparison(20E3, 1E3)
Elapsed time is 3.153404 seconds.
Elapsed time is 13.200649 seconds.
|y2-y1| = 1.5627e-07

The function is:

function Comparison(n, M)
    u = rand(n, 1);
    L = rand(M);
    gamma = rand(M, 1);

    tic
        y1 = vectorized(u, L, gamma);
    toc

    tic
        y2 = looped(u, L, gamma);
    toc

    disp(['|y2-y1| = ', num2str(norm(y2 - y1, 1))])
end

function y = vectorized(u, l, gamma)
global a Column
    M = length(gamma);
    Column = l' * gamma;

    x = bsxfun(@plus, -(1:M)', (1:length(u)) + 1);
    x(x < 1) = 1;
    a = u(x);
    a(1:M, 1:M) = a(1:M, 1:M) .* triu(ones(M));
    a = a';
    m = 1:size(a,1);
    y = arrayfun(@VectorY , m)';
end

function yi = VectorY(j)
global a Column
    yi = a(j,:) * Column;
end

function y = looped(U, l, gamma)
    n = length(U);
    M = length(gamma);

    u = U';
    L = l';
    y = zeros(n, 1);

    for i=1:n
        y(i) = [u(i:-1:max(1,i-M+1)), zeros(1,(-i+M)*(i-M<0))] * L * gamma;
    end
end

how does your vectorized solution scale in terms of memory consumption? — m.s., May 31 '15 at 14:55
This is the source of the new information: http://stackoverflow.com/questions/12522888/arrayfun-can-be-significantly-slower-than-an-explicit-loop-in-matlab-why — user1543042, May 31 '15 at 18:57

Vectorize MATLAB for loop

2 Answers2

Discussion and solution code

Benchmarking