Marginalize a table over given parameters using OpenCL

Question

Here it was asked how to generalize the procedure of marginalizing a table over given dimensions. The solution given here works very well.

I would like to have a parallelized version of that solution. In particular, I was wondering if it would be possible to use the GPU via OpenCL (I have a AMD Radeon Pro Vega 56).

In order to reduce the problem, perhaps one could initially consider the following simplified task. Let us define a test table tab:

{n1, n2, n3} = {10, 20, 30};
tab = Table[{i1, i2, i3, i1 + i2 + i3}, {i1, n1}, {i2, n2}, {i3, n3}];

Then, this command marginalizes over the first and second dimensions:

newtab = tab[[1, 1, All, {3, 4}]];
Do[newtab[[i3, 2]] = -2 Log[Total[Exp[-tab[[All, All, i3, 4]]/2],Infinity]];, {i3, n3}];

while the following one marginalizes over the second dimension only:

newtab = tab[[All, 1, All, {1, 3, 4}]];
Do[newtab[[i1, i3, 3]] = -2 Log[
  Total[Exp[-tab[[i1, All, i3, 4]]/2], Infinity]];, {i1, n1}, {i3,
n3}];

How to efficiently perform the latter operations with OpenCL?

Answer 1

Before resorting to the nuclear option one can obtain 100-fold speedup with conventional methods. More precisely, we can use Compile to guarantee that we get a PackedArray of machine precision numbers as tab2 (it also speeds up the creation of the array) and operate on this array only by vectorized methods.

{n1, n2, n3} = {10, 20, 30} 10;

tab = Table[{i1, i2, i3, i1 + i2 + i3}, {i1, n1}, {i2, n2}, {i3, n3}]; 
newtab = tab[[All, 1, All, {1, 3, 4}]];
Do[newtab[[i1, i3, 3]] = -2 Log[ Total[Exp[-tab[[i1, All, i3, 4]]/2], Infinity]];
  , {i1, n1}, {i3, n3}]; // AbsoluteTiming // First

tab2 = Compile[{{n1, _Real}, {n2, _Real}, {n3, _Real}},
     Table[{N@i1, N@i2, N@i3, N[i1 + i2 + i3]}, {i1, 1., n1, 1.}, {i2,
        1., n2, 1.}, {i3, 1., n3, 1.}]
     ][n1, n2, n3];

newtab2 = tab2[[All, 1, All, {1, 3, 4}]];
newtab2[[All, All,3]] = -2. Log[
     Total[
       Exp[-0.5 tab2[[All, All, All, 4]]], 
       {2}
       ]
     ]; // AbsoluteTiming // First

newtab == newtab2

12.321

0.110804

True

Moreover, the success of porting this to a GPU (success in the sense of a performance gain) is sensitive to wether you are planning to generate the input data (e.g., tab) on the GPU (and wether you have to get the output back from there, but that's less an issue in this case since the output is much smaller). In particular, transferring input of this size to the GPU for doing so few number crunching will hardly pay off. The following illustrates how slow the communication between CPU and GPU can be:

Needs["OpenCLLink`"]
{n1, n2, n3} = {10, 20, 30} 10;
tab3 = Compile[{{n1, _Integer}, {n2, _Integer}, {n3, _Integer}}, 
    Table[{i1, i2, i3, i1 + i2 + i3}, {i1, 1, n1}, {i2, 1, n2}, {i3, 1, n3}]
    ][n1, n2, n3];

cltab = OpenCLMemoryLoad[tab3, Integer, "Platform" -> 1, "Device" -> 3]; // AbsoluteTiming // First
tab3 = OpenCLMemoryGet[cltab]; // AbsoluteTiming // First
OpenCLInformation[$OpenCLPlatform, $OpenCLDevice, "Name"]

0.107119

0.107834

"AMD Radeon R9 M370X Compute Engine"

So, solely sending the data to the GPU is as slow as the actual computation on the CPU.