How to apply a function to each column of a 3D array?

I have a function that works on a vector (call it "test_f") and gives as output a vector. I want to apply it to each column of a 3D array ( example A= randn(10,10,10)), without using a loop. Is it possible?

[Massimo]

This is how it works with a for loop. I would like to do it without for loop, in order to implement later in a GPU.

 A=randn(10,10,10);
 for k=1:10
     for l=1:10
         B(:,k,l)=test_f(A(:,k,l));
     end
 end

Answer 1

[Joss Knight]

There isn't anything supported for gpuArray that can take any generic user function in this way. If test_f contains operations supported by pagefun then you could break it down into multiple calls to that. Or convert your array to a cell array ( mat2cell(A, 10, ones(10,1), ones(10,1) ) and process it using cellfun.

[Massimo]

To make it simple, suppose that the function test_f simply gives as output a vector with the first three elements of the input vector

function vout=test_f(vin)
vout=vin(1:3);

I would like to apply it column-wise to the array A=randn(10,10,10) that would give at the end an array with size (3x10x10).

1) You propose using mat2cell, then

B=mat2cell(A, 10, ones(10,1), ones(10,1); % 1x10x10 cell
C=cellfun(@test_f,B,'UniformOutput',false); % 1x10x10 cell
D=cell2mat(C)    %ERROR CELL2MAT does not support cell arrays containing cell arrays or objects

how can I obtain finally the 3x10x10 array?

2) You propose using pagefun. What is the right syntax of pagefun in this case?

[Joss Knight]

You wouldn't use pagefun for this, you would use indexing:

D = A(1:3,:,:);

so it's not a very good example. If your function was, say, the vector dot product between the first 5 and second 5 values in each column, you'd do that using basic array arithmetic:

D = sum( A(1:5,:,:) .* A(6:10,:,:) );

Many algorithms can use this approach. pagefun would come in if you wanted to do some matrix arithmetic. Let's imagine you want to use mtimes to compute the outer product between the columns of the array, perhaps to compute the autocorrelation. You'd reshape them appropriately and call pagefun:

A1 = reshape(A, 10, 1, 10, 10);
A2 = reshape(A, 1, 10, 10, 10);
D = pagefun(@mtimes, A1, A2);

Of course, you wouldn't actually use mtimes to do this because automatic scalar expansion does this.

D = A1 .* A2;

If you have a function that is too complicated to be broken down like this. You can try the cellfun approach.

Acell =  mat2cell(A, 10, ones(10,1), ones(10,1));
D = reshape( cellfun(@test_f, Acell), [], 10, 10 );

But there is no parallelism here, so the array passed to each function is going to have be large to make it worth using the GPU. Since the above is no different to a loop, you're back where you started.

[Massimo]

So there is no advantage to implement in this case the vectorization in the GPU for a function that works on vectors. In fact using a for loop the calculation is much faster than using mat2cell in the CPU and faster than using mat2cell in the GPU. As an example:

tic
C=randn(100,100,100);
M_OUT=ones(3,100,100);
for x=1:100
    for y=1:100
        M_OUT(:,x,y)=test_mat(C(:,x,y));
    end
end
toc   %Elapsed time is 0.059773 second

and

tic
A=randn(100,100,100);
B= mat2cell(A, 100, ones(100,1), ones(100,1))  ;
C=cellfun(@test_mat,B,'UniformOutput',false);
D=[C{:}];
Z=reshape(D, [3,100,100]);
toc % Elapsed time is 0.321378 seconds.

and

tic
A=randn(100,100,100,'gpuArray');
B= mat2cell(A, 100, ones(100,1), ones(100,1))  ;
C=cellfun(@test_mat,B,'UniformOutput',false);
D=[C{:}];
Z=reshape(D, [3,100,100]);
toc %Elapsed time is 1.884156 seconds.

[Joss Knight]

Your expectations for the capabilities of a GPU are misguided in this case. There's almost nothing a GPU can do on 100 values faster than the CPU. You need to give it more data. One way is to vectorize your code, which means working out how to formulate the equations so that all the data is processed at once. I can't help you do that as long as I've no idea what test_mat is doing.