# Marginalize a table over given parameters using OpenCL

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?

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.

• thanks, usually I do use PackedArray's, I put the command ToPackedArray before feeding the tab` to the function. Your solution is anyway more elegant (and 4 times faster for this example)! Regarding the GPU, the table is usually obtained after solving a bunch of stuff (including differential equations) and so I believe it would be complicated to do it directly on the GPU. I understand the lesson. What are the most efficient uses of a good GPU? Apr 23, 2018 at 20:25
• Well, I am not an expert in GPU programming. There are definitely valuable applications, e.g. in image processing where each process needs access only to a small part in an image. People use GPUs also a lot in machine learning but the price they have to pay is that they have to use lousy algorithms for solving the optimization problems. Apr 23, 2018 at 20:39
• But it' s not only the communication between CPU and GPU that makes life hard. For example, one has to make sure that the global address space on the GPU (a huge vector) is read in accending order; otherwise you wil experience severe slowdown. Random access is one of the many things that a GPU can' t do well. That makes GPUs nearly useless for PDEs on unstructured grids and that' s also the reason why I did not dive deeper into this topic. Apr 23, 2018 at 20:39