title: Why is implicit expansion not always faster on GPU?
body: <p>I am trying to subtract two matrices of "different" sizes in Matlab on my GPU. </p>
  
  <pre><code>a = randn(87,1,10,'single','gpuArray');
  b = randn(87,100000,1,'single','gpuArray');
  
  dev = gpuDevice;
  tic
  out1 = a-b;
  wait(dev)
  toc
  
  tic
  a = repmat(a,[1 size(b,2) 1]);
  b = repmat(b,[1 1 size(a,3)]);
  out2 = a-b;
  wait(dev)
  toc
  </code></pre>
  
  <p>As you can see, there is implicit expansion along the second and third dimensions. Typically, implicit expansion is much faster than replicating matrices. However, in this case it gives the result:</p>
  
  <pre><code>Elapsed time is 0.033324 seconds.
  Elapsed time is 0.018052 seconds.
  </code></pre>
  
  <p>This is perplexing because replicating matrices is much faster than just letting implicit expansion handle the difference. I tested this on the CPU and got: </p>
  
  <pre><code>Elapsed time is 0.078726 seconds.
  Elapsed time is 0.356435 seconds.
  </code></pre>
  
  <p>What exactly is going on here? Why is implicit expansion so much slower?</p>
  
  <p>I tested this on Matlab 2018b and 2019a prerelease and got this result on both.</p>
  
  <p>Edit: Due to popular demand I used the timeit and gputimeit functions on both and this just reaffirmed my results</p>
  
  <pre><code>clear
  close all
  a = randn(87,1,10,'single','gpuArray');
  b = randn(87,100000,1,'single','gpuArray');
  
  dev = gpuDevice;
  disp(gputimeit(@() a-b))
  
  disp(gputimeit(@() a&gt;b))
  
  t1 = gputimeit(@() repmat(a,[1 size(b,2) 1]))+gputimeit(@() repmat(b,[1 1 size(a,3)]));
  a = repmat(a,[1 size(b,2) 1]);
  b = repmat(b,[1 1 size(a,3)]);
  
  disp(t1+gputimeit(@() a-b))
  
  disp(gputimeit(@() a&gt;b))
  
  clear
  close all
  a = randn(87,1,10,'single');
  b = randn(87,100000,1,'single');
  
  disp(timeit(@() a-b))
  
  disp(timeit(@() a&gt;b))
  
  t1 = timeit(@() repmat(a,[1 size(b,2) 1]))+timeit(@() repmat(b,[1 1 size(a,3)]));
  a = repmat(a,[1 size(b,2) 1]);
  b = repmat(b,[1 1 size(a,3)]);
  
  disp(t1+timeit(@() a-b))
  
  disp(timeit(@() a&gt;b))
  </code></pre>
  
  <p>And the times are: </p>
  
  <pre><code>    0.0249880485773384
  
      0.0191386963821538
  
      0.0159509649972431
  
     0.00518230250048866
  
      0.0980248506300088
  
      0.0545817410501994
  
       0.283484452903233
  
      0.0540456719899175
  </code></pre>