Filling in a 2-D matrix slice by slice with parfor

Hi,

I have a 2-d matrix where one dimension is a unique starting state, the second is time, and I need to compute entries in this matrix. Given the starting state, the evolution in time follows a recursive formula. For instance:

A = zeros(1000,1000);
A(:,1) = 1;
for i = 1:1000
     for t = 2:1000
          A(i,t) = 1 + A(i,t-1);
     end
end

In the actual application there is no way around the recursive formulation in the t dimension. I would like to replace the outer loop with parfor, but I get an error since A(i,t-1) appears on the RHS of the assignment operator. Conceptually the slices in the i dimension are independent so it should be possible to parallelize.

Any advice on how to implement parfor?

Answer 1

[Matt J]

It's hard to imagine PARFOR having any benefit on data of such a small size. Can't you just use filter()? It's already vectorized column-wise and, as a .bi function, probably already has lots of internal parallelization.

[Daniel]

In my code the inner loop takes about 2.5 seconds to execute due to some big kron operations. The outer loop has about 2000 iterations. It would nice if it were only a few minutes instead of >1 hour to run

I've never used filter() before and the help file notation is a bit general. Could you show me what you mean on the example above?

[Matt J]

Your example could be implemented using FILTER as follows

 >> X=ones(10); A=filter(1,[1,-1],X).'

It could also be implemented with only one for-loop as follows

     for t = 2:1000
          A(:,t) = 1 + A(:,t-1);
     end

I don't know if any of this helps anymore because I'm not sure if your recursion is linear or if the recursion rule is i-dependent. You might need to refine your example, showing more fully what its complications are.

[Daniel]

It's i-dependent in that it'll need to access i-dependent values from another fixed matrix, and it's definitely linear.

The actual also operates on a 3 dimensional matrix, with a vectorization over the 3rd dimension. Something like

A = zeros(2000,100,50);
eps = randn(2000,100,50);
A(:,1,1) = 1;
for i = 1:1000
     for t = 2:100
          A(i,t,:) = 1 + 0.5*A(i,t-1,:) + eps(i,t,:);
     end
end

[Matt J]

I don't see any kron operations in there. All your examples so far are within the scope of FILTER. The third dimension doesn't matter. You can use RESHAPE to convert to a 2D array and back.