How to left divide 4D arrays

If I have two 2D gpuArrays f(NxM) and f1(NxM) of same size, then I can left divide them like this:

beta=arrayfun(@(n) f(:,n)\f1(:,n), 1:size(f,2));

What would be the way to do this if both f and f1 are 4D arrays of size (NxMxZxK)?

Answer 1

[Walter Roberson]

beta = arrayfun(@(n) f(:,n)\f1(:,n), 1:numel(f)/size(f,1));

after which you will need to reshape() to the appropriate size.

[Mantas Vaitonis]

Great, thank You.

[Mantas Vaitonis]

@ Walter Roberson, maybe there is way to speed up this code? Because now when I implement it to my algorithm, it become the bottleneck and takes 95% of all algorithm. I use gpuArrays, and maybe delay is because it saves result to cell arrays, which are slower on gpu?

[Walter Roberson]

That code does not save the result to cell arrays. Each call to arrayfun with those arguments would result in a 1 x 1 scalar numeric value.

[Mantas Vaitonis]

When I implement this code it results in cell array 1xK with cells size 1x1, what about the possibility to speed up?

[Mantas Vaitonis]

I use f and f1, both are 4D gpuArrays, then I use this code:

 beta2= arrayfun(@(n) f(:,n)\f1(:,n), 1:numel(f)/size(f,1),'uni',0);

Result beta2 is cell array. If I change numel(f)/size(f,1) from double to gpu double:

fd=gpuArray(numel(f)/size(f,1));

 beta2= arrayfun(@(n) f(:,n)\f1(:,n), 1:fd);

I get this error: Function passed as first input argument contains unsupported or unknown function 'colon'. Which is true, because arrayfun does not work with indexing.

[Walter Roberson]

Notice that my version of the code does not have 'uni', 0.

>> clearvars;aaa = randi([-10 10], 3,4,5,6);
>> bbb = randi([-10 10], 3,4,5,6);
>> f = gpuArray(aaa);
>> f1 = gpuArray(bbb);
>> beta = arrayfun(@(n) f(:,n)\f1(:,n), 1:numel(f)/size(f,1));
>> size(beta)
ans =
     1   120

No 'uni', 0 needed.

[Mantas Vaitonis]

@ Walter Roberson if I use your code without the 'uni', 0 I get this error "Error using arrayfun gpuArray output type is not supported. Set 'UniformOutput' to false." my version of Matlab R2016b.

[Walter Roberson]

Ah, it appears to be R2017b that added the needed support.

Perhaps you could get somewhere by reshaping the array to 3D and using pagefun()

[Mantas Vaitonis]

I know that I am tempting with your patience, but what would be the option with pagefun, because as I know it runs only with build in functions and I did try rdivide/ldivide, I can reshape even to 2D by using your suggestion:

f2d=f(:,numel(f)/size(f,1));

And what would be further steps?

[Walter Roberson]

The / operator is formally named mrdivide. The R2016b pagefun documentation lists @mrdivide specifically as being supported.

[Mantas Vaitonis]

I did try your code on R2017b and yes it does work without 'uni' 0, however it did not change the speed of calculations, if you try with bigger 4D arrays it gets slow very fast:

clear;
aaa = randi([-10 10], 200,4,10,100 ); 
bbb = randi([-10 10], 200,4,10,100 ); 
f = gpuArray(aaa ); 
f1 = gpuArray(bbb ); 
tic
beta = arrayfun(@(n) f(:,n)\f1(:,n), 1:numel(f)/size(f,1 )); 
toc

Is this the only way possible?

[Walter Roberson]

I do not know what rate you expect for those. That takes 4 to 6 seconds on my system:

Elapsed time is 5.437293 seconds.
Elapsed time is 3.842687 seconds.
Elapsed time is 3.824775 seconds.

[Mantas Vaitonis]

On my machine it takes about 7,5 seconds:

Elapsed time is 7.577014 seconds.

I use other functions with arrayfun, pagefun on same algorithm with same data and it for evey other function it takes less then 0,01 second. One more thing when I use profiler it says there were 4000 calls in der>@(n)f(:,n)\f1(:,n), is in tuppose to pas data in parallel? I have quite a big chunks of data and looking for ways to speed the calculations.

[Walter Roberson]

The left division of two real-valued column vectors A\B, is equal to the dot product of the two vectors divided by the square of the norm of A (that is, the sum of squares of the elements of A)

[Mantas Vaitonis]

@Walter Roberso, great thank you, using the method you did suggest speed up is more then 100x times:

clear;
aaa = randi([-10 10], 5,4,10,100 ); 
bbb = randi([-10 10], 5,4,10,100 ); 
f = gpuArray(aaa ); 
f1 = gpuArray(bbb ); 
tic
ff=f;
ff1=f1;
r=dot(ff,ff1);
rr=ff.^2;
S = sum(rr,1);

 ss=arrayfun(@rdivide,r,S);