Many small Eigenvalue Decompositions in parallel on GPU?

I have some code that involves a couple billion 3x3 and 4x4 eigenvalue decompositions. I have run this code with parfors on the CPU and the runtime is just barely bearable, but I'd really like to speed this up.

I have a GTX 780 available. I realize that a GPU is generally better suited for large matrix operations than a large number of small matrix operations. I looked at pagefun, which looks like the best way that Matlab has to run many small matrix operations in parallel. However, the functions available for pagefun are all element by element operations, with a few exceptions such as mtimes, rdivide, and ldivide. Unfortunately eig is not one of those functions.

Is there any other way to run this code on the GPU?

[Matt J]

Are you sure you mean "several thousand"? My old machine from 2008 can do 10000 such decompositions without breaking a sweat,

>> tic; for i=1:10000, eig(rand(4)); end; toc
Elapsed time is 0.196188 seconds.

[ervinshiznit]

Oops. I just said "several thousand" without actually looking at how many times I'm calling eig. Looking at it, it's actually 2,200,570,000 calls to eig.

I'll edit the original post

Of course this code involves other calculations as well which contribute to the runtime, but the eig is the slowest portion.

Answer 1

[Brian Neiswander]

The "pagefun" function does not currently support the function "eig". However, note that the "eig" function will accept GPU arrays generated with the "gpuArray" function:

X = rand(1e3,1e3);
G = gpuArray(X);
Y = eig(G);

Depending on your data, this can be faster than the non-GPU approach but it is not parallelized across the pages.

It is possible to implement your own CUDA kernel using the CUDAKernel object or MEX functions. This allows for you to process custom functions using a distribution scheme of your choice. See the links below for more information:

http://www.mathworks.com/help/distcomp/parallel.gpu.cudakernel.html

http://www.mathworks.com/help/distcomp/run-cuda-or-ptx-code-on-gpu.html#bsic5ih-1

[ervinshiznit]

I already tried gpuArray. It's far too slow, the transfer times to and from the GPU kill me. It does provide a speedup for larger matrices, but not 3x3 or 4x4.

CUDA kernels will not work for me because that's a lot of development time that I do not have. Looks like I'm just stuck with the runtimes.

[Birk Andreas]

So, its already 2019 and there are already some MAGMA eigenvalue functions implemented. However, still no eig for pagefun...

What prevents the progress?

Could you give an estimate, when it will be implemented?

It would really be very welcome!

Answer 2

[Joss Knight]

Have you tried just concatenating your matrices in block-diagonal form and calling eig? You may then be limited by memory, but the eigenvalues and vectors of a block-diagonal system are just the union of the eigenvalues and vectors of the blocks:

N = 1000;
A = rand(3,3,N);
maskCell = mat2cell(ones(3,3,N),3,3,ones(N,1));
mask = logical(blkdiag(maskCell{:}));
Ablk = gpuArray.zeros(3*[N,N]);
Ablk(mask) = A(:);
[Vblk,Dblk] = eig(gpuArray(Ablk));
V = reshape(Vblk(mask), [3 3 N]);
D = reshape(Dblk(mask), [3 3 N]);

You should then find that A(:,:,i)*V(:,:,i) == V(:,:,i)*D(:,:,i) as required. Because of the way eigendecomposition works, I would expect the extra unnecessary zeros not to affect the performance much, the system should converge straightforwardly and parallelize well.

[ervinshiznit]

This is a great idea. I'll test this out.

Furthermore, this would let me scale to matrices greater than 4x4, which is something I'd like to be able to do.

[ervinshiznit]

I tried it out. Memory issues aside (I could break this up into chunks to achieve it), sadly this is still slower on the GPU than the CPU.

[Joss Knight]

What is the bottleneck? Is eig itself slower on the GPU than the CPU? Run

Gblk = gpuArray(Ablk);
timeit(@()eig(Ablk),2)
gputimeit(@()eig(Gblk),2)

[Joss Knight]

Also, I see that the GTX 780 has a terrible double-precision performance of 166 GFlops versus 3977 for single precision. Try running your code in single precision.

[kunx]

thank you. your direction is very helpful.

Answer 3

[James Tursa]

If you just need the eigenvalues, you might look at this FEX submission by Bruno Luong:

http://www.mathworks.com/matlabcentral/fileexchange/27680-multiple-eigen-values-for-2x2-and-3x3-matrices

Maybe you can expand it for 4x4 as well.

[ervinshiznit]

I also need the eigenvalues.

You gave me an idea though - what if I explicitly wrote out the formula for solving a 3x3 system (not as a matrix inversion, but 3 equations for solving each coordinate of a 3x1 vector individually in the 3x3 case).

That's not something that easily lends itself to a GPU implementation, but it may still be faster to run that on the CPU than calling eig for small enough matrices (like my 3x3's and 4x4's).

[Joss Knight]

You don't need to use an explicit formula - pagefun supports mldivide (the backslash \ operator) for solving small systems.

[ervinshiznit]

I know, but like I said in a comment to Brian's answer, transfer times of 3x3 and 4x4 matrices to the GPU kill me. I was saying that maybe I should do an explicit formula on the CPU, not the GPU. But your answer of doing a block diagonal matrix might work out.

[Joss Knight]

Why do you need to transfer 3x3 and 4x4 matrices to the GPU independently? Just transfer it all as one 3D array. You have to anyway to use pagefun.