How do I use @mldivide together with pagefun on the GPU?

I need to solve a lot of batches of xA = B problems and I'm experimenting with different solutions to see which type is fastest. And since pagefun supports mldivide I suspect it could be the perfect function for what I want to do. However, I have never fully understood the concept behind how it works and the documentation only gives examples of multiplication problems that I can't completely relate to what I want to do.

So could someone please just provide a small example or tell me how to use pagefun with mldivide to solve several xA = B at the same time so I can see how it should be set up?

Say for example I have an A of size = [100 10 15] and a B of size = [100 1 15] and I want to solve the problem for the 15 different versions along the third dimension. So in a loop it would be:

for i = 1 : 15
A(:,:,i)\B(:,1,i) 
end

But how do I use pagefun to do this in parallel without the loop? If I write it as:

A = gpuArray.rand(100, 10,15);  
B = gpuArray.rand(100, 1, 15);
x = pagefun(@mldivide, A, B);

I get:

Error using gpuArray/pagefun
The number of rows and columns of each page of the divisor must be equal when
using MLDIVIDE or MRDIVIDE with PAGEFUN.

This is probably very trivial to everyone else, but I just don’t get it. Could someone clarify how this is done? Thanks!

Answer 1

[Joss Knight]

Your systems are rectangular. You can only solve square systems with pagefun, sorry.

[Peta]

Damn. Do you know if there is an alternative function that can solve this in parallel on the GPU then that I should try to use instead? I’ve tried placing my A’s in a structure and then use structfun, which works but the performance isn’t that amazing.

[Joss Knight]

You just want a pseudo inverse like inv(A'*A)*A'. So why not try

At = pagefun(@ctranspose, A);
AtA = pagefun(@mtimes, At, A);
pseudoInv = pagefun(@mldivide, AtA, At);
X = pagefun(@mtimes, pseudoInv, B);

Not as robust as a proper QR-based solve but could be good enough for you.