bugGNU Octave - Bugs: bug #59840, repmat and repelem slower than...

 
 

bug #59840: repmat and repelem slower than needed, faster implimentation suggested

Submitter:  None
Submitted:  Fri 08 Jan 2021 01:49:29 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Performance
Status:  In Progress Assigned to:  None
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  * Open Release:  * 6.1.0
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 30 Jan 2023 02:42:36 PM UTC, comment #11: 

I haven't seen significant changes to repmat or repelem in some time.  Note that when I helped write Octave's repelem a few years ago, the main focus was creating a compatible implementation, catching all of the corner cases.  I did some small optimizations around small inputs, and avoided adding the overhead of repmat calls for obvious reasons, but didn't do extensive speed tests. There was, and still is, an existing Octave-only function called repelems that only does 1D repeats and only works on numeric arrays, but as a compiled function it's much faster. A lot of Octave code still uses it over repelem because of the speed difference.

The extensive look into problem-size algorithm choice is generally useful.  We ran into a similar discussion trying to decide on the appropriate algorithm for 2D triangulation (see delaunayn bug #60818).  It turns out one algorithm was faster up to complexity X, another kron-based one was better beyond that but resulted in much higher memory use to the point where a number of places it was faster couldn't be completed due to OOM issues. Some lightweight algorithm selection routine could be generally useful. But memory footprint of the approach really should be part of that consideration as well and would require machine specific knowledge.  (maybe OOM is less of an issue here, though since the output is intended to be a multiple in size of the input, if any intermediates don't exceed this significantly)

Nicholas Jankowski <nrjank>
Group Member
Mon 30 Jan 2023 04:20:09 AM UTC, comment #10: 

What happened to this? Is it adopted?

Anonymous
Tue 19 Jan 2021 05:13:54 PM UTC, comment #9: 

Given that kron requires the multiplication of each A(i,j) by the matrix B it had been my thought that the product of dimensions would be the right criteria.  That's why I was surprised to find there was a dependence on the actual values of m, n, p, q rather than just their product.  I will review your results on my computer when I have the chance.

One thing to note is that there may be much more variability in the kron-based approach compared to the indexing approach.  The indexing approach relies solely on creating large contiguous blocks of memory (essentially malloc()) and then copying data (essentially memcpy).  There will be a dependence on processor bus speeds and the speed of attached memory, but that's about it.

The kron-based approach uses multiplication which, for starters, may be software-emulated on very simple processors versus having a dedicated hardware multiplier.  Even when a HW multiplier is present, there may be varying number per processor.  This may mean it is too difficult to find a cut-off that works well for a range of machines.  But, I will take a look.

Rik <rik5>
Group administrator
Tue 19 Jan 2021 04:36:46 AM UTC, comment #8: 

After letting it run all night and generating over a gigabyte of timing for every combination data I found a simple function for where to switch between the kron and idx implementations for both repmat and repelem:

if m*n*p*q<128e3 and double precision, kron normally faster otherwise idx faster

if m*n*p*q<128e3*2 and single precision, kron normally faster otherwise idx faster

There is a significant increase in computation time at the above switch location for both methods.

As having the precision doubled the product constant I'm confident the difference in speed is related to memory. This may mean the constant is different for different computers.

I have attached the code and some small data files. I have drawn a line for the transition between the methods below the timing plots, currently it is for single precision (remove the times 2 for double precision).

(file #50748)

Anonymous
Mon 18 Jan 2021 02:11:53 PM UTC, comment #7: 

To look for a better function to switch implementations I decided to plot the four dimensional domain to get an idea where one function was better than the other. I only found one case where the current implementation of repmat was clearly better than the kron implementation.

let m and n be the matrix size, p the number of times vertically repeated and q be the number of times horizontally repeated.

The current implementation is better when p=1, q>12 m>20, n>20 (very rough) also as q increases m and n decrease from the above numbers, otherwise kron is generally faster by a tiny amount. I plotted 10 points between 1 and 100 in each variable (10000 points total) and ran each repmap case 10 times. Total run time to create the data was around two hours and overall the kron way was around 24 seconds faster.

I have not run repelem yet.

I have attached my data (repmattimekronindx), the method to visualize and create the data (repmatrepelemfaster.m). To create new data change the first line to true.

(file #50743, file #50744)

Anonymous
Sun 17 Jan 2021 02:05:53 AM UTC, comment #6: 

I wrote an objective function to attempt to find a constant to use to switch between an indexing implementation or a kron-based implementation.  The function is attached and listed below.


function retval = objf (x)

  m = round (x(2));  n = round (x(1));

  p = 50;  q = 50;
  A = rand (m, n);
  N = 1000;

  tic;
  for i = 1:N
    tmp = repmatidx (A, p, q);
  endfor
  tm.idx = toc;

  tic;
  for i = 1:N
    tmp = repmatkron (A, p, q);
  endfor
  tm.kron = toc;

  retval = tm.idx - tm.kron;

endfunction


First, I found that individual measurements of timing bounced around too much so I ran each timing test 30 times (N = 30) in order to try to get more consistent results.  That still didn't work very well.  Eventually I bumped it to N = 1000.

Using fminsearch() did not work for me.  The algorithm had difficulty converging.  This may be because there is still to much  variability in the result returned from the objective function.

But I did notice some trends.  When M and N, and likewise P and Q, are roughly equal implying roughly square matrices the indexing operations tend to be faster.  If, instead, M and N or P and Q are vastly different implying a very rectangular matrix then kron works better.

The choice to switch between the two implementations probably needs to be more in depth than just summing the number of rows and columns.

As an example, with p=25, q=25:

objf ([100, 100]) = -1.4   # implies indexing is better
objf ([199, 1  ]) = +.06   # implies kron is better

In both of these case M+N = 200 so a decision based only on that sum cannot distinguish these two cases.

It does make me think that this is very input data dependent.  Can you discover a more precise criterion for switching between implementations?

(file #50727)

Rik <rik5>
Group administrator
Sun 10 Jan 2021 04:29:29 AM UTC, comment #5: 

Just noticed I forgot to put the loop counter back to 1000 for the loop solution. Now it's the fastest to about 30 iterations, then it's the slowest by a huge amount.

Anonymous
Sun 10 Jan 2021 01:31:30 AM UTC, comment #4: 

After a significant amount of additional time on this, I have addressed each comment below and provided a new implementation that is generally faster than the existing method with all of the features of the existing method. I found for large inputs the existing method is faster than the kron method but the kron method is faster for small inputs. For both repmat, repelem I copied the existing code and modified the relevant section. I put a row of comment signs above and below the section I modified. The modifications are a few changes to the checks in repmat and repelem with the inclusion of kron into repelem.

The loop solution is still the fastest for my particular problem but appears to be slower if I put it into a function so I think the speedup is related to avoiding copying memory in and out of a function. It is still faster for a few hundred iterations and as the rows were variables for repeated single variable integration I will never have large numbers of rows. The reason for the difference in output in my original submission was a mistake introduced after my copy and paste, there should have been a +1 after size on the first line. I clearly deleted it after copying my original code as it was present where I copied it from and I did not check the output before submitting the bug.

I worked out a simple function for were to switch between the kron method and the existing method by creating an objective function with inputs of the size of the matrix to be repeated and the number of times to be repeated in each direction. The objective function then evaluated the time using both methods and returned the absolute difference in time. I minimized this function using fminsearch (natural number inputs so the function is not smooth) with a uniform random guess between 1 and 100 for each variable, 100 times (if you feel like spending more time you can do more). Using this output I decided if the sum of the inputs (matrix size plus the number of times repeated in each direction) is less than 175 use the kron method otherwise use the existing method. I chose this as the mean sum of inputs where the evaluation time was the same was around 175, the standard deviation was around 50.

My code is below: it contains the optimization code which takes a few minutes to run, the original tests which can now easily be extended to more cases, and the proposed replacement functions repmatnew and repelemnew. repmatcmpspeed and repelemcmpspeed are identical excluding the sum of input checks which is set so the kron method always runs to perform the optimization.

Another thought with sparse matrixes and repelem perhaps the kron method should always be used to keep it consistent with repmat. I have not checked the speed difference in this case.

I have also attached plots of the input sizes where the new and old functions have the same evaluation time so the slow first section of the code can be skipped.


function repmatrepelemfaster

a=(1:100);
xx=nan(4,length(a));
for k=a
  k
  [x,r,ex]=fminsearch(@repmatswitchmeth,100*rand(4,1));
  if ex==1
    xx(:,k)=abs(round(x));
  end
end
close all;
figure;
plot(a,xx(1,:),'o',a,xx(2,:),'s',a,xx(3,:),'p',a,xx(4,:),'d',a,sum(xx,1));
legend('m','n','a','b','sum')
title('repmat');
inputsum=sum(xx,1);
inputsum(isnan(inputsum))=[];
mean(inputsum)
std(inputsum)

a=(1:100);
xx=nan(4,length(a));
for k=a
  k
  [x,r,ex]=fminsearch(@repelemswitchmeth,100*rand(4,1));
  if ex==1
    xx(:,k)=abs(round(x));
  end
end
figure;
plot(a,xx(1,:),'o',a,xx(2,:),'s',a,xx(3,:),'p',a,xx(4,:),'d',a,sum(xx,1));
legend('m','n','a','b','sum')
title('repelem');
inputsum=sum(xx,1);
inputsum(isnan(inputsum))=[];
mean(inputsum)
std(inputsum)







q=reshape(1:6,[],2)

repmat(q,3,1)
repmatnew(q,3,1)

repelem(q,3,2)
repelemnew(q,3,2)

ab=[1 50;
50 1;
50 50];
qq=cell(3,1);
qq{1}=1:10;
qq{2}=qq{1}';
qq{3}=reshape(1:100,10,[]);

for i=1:length(qq)
for j=1:size(ab,1)
  a=ab(j,1);
  b=ab(j,2);
  q=qq{i};
  fprintf('b=%d, c=%d, size q=%d,%d\n',a,b,size(q))
  disp('inbuilt repmat');
  tic
  for k=1:1000
    repmat(q,a,b);
  end
  toc
  disp('repmat using kron')
  tic
  for k=1:1000
    repmatnew(q,a,b);
  end
  toc
  disp('inbuilt repelem')
  tic
  for k=1:1000
    repelem(q,a,b);
  end
  toc
  disp('repelem using kron')
  tic
  for k=1:1000
    repelemnew(q,a,b);
  end
  toc
end
end
fprintf('What I was doing when I observed the issue:\nInsert a row into a column vector repeting the elements of the vector\nThis was used to plot a slice through a multivariable function.\nSpeed became an issue when I tried to intergrate the function.\n');

np=150;
xi=rand(1,np);
xo=rand(20,1);
i=2;
disp('Origional attempt')
tic
for l=1:10000
  x=[repmat(xo(1:i-1),1,np);xi;repmat(xo(i:end),1,np)];
end
toc
disp('with faster repmat')
tic
for l=1:10000
  x=[repmatnew(xo(1:i-1),1,np);xi;repmatnew(xo(i:end),1,np)];
end
toc
disp('single repmat call')
tic
for l=1:10000
  x=zeros(size(xo,1)+1,size(xi,2));
  x([1:i-1 i+1:end],:)=repmat(xo,1,np);
  x(i,:)=xi;
end
toc
disp('with faster repmat')
tic
for l=1:10000
  x=zeros(size(xo,1)+1,size(xi,2));
  x([1:i-1 i+1:end],:)=repmatnew(xo,1,np);
  x(i,:)=xi;
end
toc
disp('loop solution')
tic
for l=1:1
  nv=size(xo,1)+1;
  x=zeros(size(xo,1)+1,size(xi,2));
  x(i,:)=xi;
  for k=[1:i-1]
    x(k,:)=xo(k);
  end
  for k=i+1:nv
    x(k,:)=xo(k-1);
  end
end
toc

end
function r=repmatswitchmeth(x)
q=rand(max(abs(round(x(1))),1),max(abs(round(x(2))),1));
a=max(abs(round(x(3))),1);
b=max(abs(round(x(4))),1);
tic;
repmat(q,a,b);
t1=toc;
tic;
repmatcmpspeed(q,a,b);
t2=toc;
r=abs(t2-t1);
end
function r=repelemswitchmeth(x)
q=rand(max(abs(round(x(1))),1),max(abs(round(x(2))),1));
a=max(abs(round(x(3))),1);
b=max(abs(round(x(4))),1);
tic;
repelem(q,a,b);
t1=toc;
tic;
repelemcmpspeed(q,a,b);
t2=toc;
r=abs(t2-t1);
end



function x = repmatcmpspeed(A, m, varargin)

  if (nargin < 2)
    print_usage ();
  endif

  if (nargin == 3)
    n = varargin{1};
    if (! isempty (m) && isempty (n))
      m = m(:).';
      n = 1;
    elseif (isempty (m) && ! isempty (n))
      m = n(:).';
      n = 1;
    elseif (isempty (m) && isempty (n))
      m = n = 1;
    else
      if (all (size (m) > 1))
        m = m(:,1);
        if (numel (m) < 3)
          n = n(end);
        else
          n = [];
        endif
      endif
      if (all (size (n) > 1))
        n = n(:,1);
      endif
      m = m(:).';
      n = n(:).';
    endif
  else
    if (nargin > 3)
      ## input check for m and varargin
      if (isscalar (m) && all (cellfun ("numel", varargin) == 1))
        m = [m varargin{:}];
        n = [];
      else
        error ("repmat: all input arguments must be scalar");
      endif
    elseif (isempty (m))
      m = n = 1;
    elseif (isscalar (m))
      n = m;
    elseif (ndims (m) > 2)
      error ("repmat: M has more than 2 dimensions");
    elseif (all (size (m) > 1))
      m = m(:,1).';
      n = [];
    else
      m = m(:).';
      n = [];
    endif
  endif
  idx = [m, n];

  if (all (idx < 0))
    error ("repmat: invalid dimensions");
  else
    idx = max (idx, 0);
  endif

  if (numel (A) == 1)
    ## optimize the scalar fill case.
    if (any (idx == 0))
      x = resize (A, idx);
    else
      x(1:prod (idx)) = A;
      x = reshape (x, idx);
    endif
  elseif (ndims (A) == 2 && length (idx) < 3)
################################################################################
    m=rows(A);
    n=columns(A);
##    Change the 1000 to 175
    if (issparse(A)||((m+n+sum(idx)<1000)&&(isnumeric(A)||islogical(A))))
      x=kron(ones(idx),A);
    else
      ## indexing is now faster, so we use it rather than kron.
      p = idx(1); q = idx(2);
      x = reshape (A, m, 1, n, 1);
      x = x(:, ones (1, p), :, ones (1, q));
      x = reshape (x, m*p, n*q);
    endif
################################################################################
  else
    aidx = size (A);
    ## ensure matching size
    idx(end+1:length (aidx)) = 1;
    aidx(end+1:length (idx)) = 1;
    ## create subscript array
    cidx = cell (2, length (aidx));
    for i = 1:length (aidx)
      cidx{1,i} = ':';
      cidx{2,i} = ones (1, idx (i));
    endfor
    aaidx = aidx;
    ## add singleton dims
    aaidx(2,:) = 1;
    A = reshape (A, aaidx(:));
    x = reshape (A (cidx{:}), idx .* aidx);
  endif

endfunction

function retval = repelemcmpspeed(x, varargin)

  if (nargin <= 1)
    print_usage ();

  elseif (nargin == 2)

    R = varargin{1};

    if (isscalar (R))

      if (! isvector (x))
        error (["repelem: %dD Array requires %d or more input " ...
                "arguments, but only %d given"], ...
               ndims (x), ndims (x) + 1, nargin);
      endif

      if (iscolumn (x))
        ## element values repeated R times in a col vector
        retval = x.'(ones (R, 1), :)(:);
      else
        ## element values repeated R times in a row vector
        retval = x(ones (R, 1), :)(:).';
      endif

    elseif (isvector (x) && isvector (R))

      ## vector x with vector repeat.
      if (numel (R) != numel (x))
        error (["repelem: R1 must either be scalar or have the same " ...
                "number of elements as the vector to be replicated"]);
      endif

      ## Basic run-length decoding in function prepareIdx returns
      ## idx2 as a row vector of element indices in the right positions.
      idx2 = prepareIdx (R);
      ## Fill with element values, direction matches element.
      retval = x(idx2);

    else # catch any arrays passed to x or varargin with nargin==2
      error (["repelem: when called with only two inputs they must be " ...
              "either scalars or vectors, not %s and %s."],
             typeinfo (x), typeinfo (R));
    endif

  elseif (nargin == 3)  # special optimized case for 2-D (matrices)

    ## Input Validation
    xsz = size (x);
    vector_r = ! (cellfun (@numel, varargin) == 1);

    ## 1. Check that all varargin are either scalars or vectors, not arrays.
    ##    isvector gives true for scalars.
    if (! (isvector (varargin{1}) && (isvector (varargin{2}))))
      error ("repelem: R1 and R2 must be scalars or vectors");

    ## 2. check that any repeat vectors have the right length.
    elseif (any (cellfun (@numel, varargin(vector_r)) != xsz(vector_r)))
      error (["repelem: R_j vectors must have the same number of elements " ...
              "as the size of dimension j of X"]);
    endif

    ## Create index arrays to pass to element.
    ## (It is no slower passing to prepareIdx
    ## than checking and doing scalars directly.)
################################################################################
## See bug 59840 for reasons for different methods
    if (isnumeric(x)||islogical(x))&&ismatrix(x)&&(issparse(x)||sum([xsz varargin{1} varargin{2}])<1000)%Change to 175, 1sd 50
      retval=kron(x,ones(varargin{1},varargin{2}));
    else
      idx1 = prepareIdx (varargin{1}, xsz(1));
      idx2 = prepareIdx (varargin{2}, xsz(2));

    ## The ":" at the end takes care of any x dimensions > 2.
      if issparse()
        retval = x(idx1, idx2);
      else
        retval = x(idx1, idx2, :);
      endif
    end
##################################################################################
  else  # (nargin > 3)

    ## Input Validation
    xsz = size (x);
    n_xdims = numel (xsz);
    vector_r = ! (cellfun (@numel, varargin) == 1);

    ## 1. Check that all repeats are scalars or vectors
    ##    (isvector gives true for scalars);
    if (! all (cellfun (@isvector, varargin(vector_r))))
      error ("repelem: R_j must all be scalars or vectors");

    ## 2. Catch any vectors thrown at trailing singletons,
    ##    which should only have scalars;
    elseif (find (vector_r, 1, "last") > n_xdims)
      error ("repelem: R_j for trailing singleton dimensions must be scalar");

    ## 3. Check that the ones that are vectors have the right length.
    elseif (any (cellfun (@numel, varargin(vector_r)) != xsz(vector_r)))
      error (["repelem: R_j vectors must have the same number of elements " ...
              "as the size of dimension j of X"]);

    endif

    n_rpts = nargin - 1;
    dims_with_vectors_and_scalars = min (n_xdims, n_rpts);

    ## Preallocate idx which will contain index array to be put into element.
    idx = cell (1, n_rpts);

    ## Use prepareIdx() to fill indices for dimensions that could be
    ## a scalar or a vector.
    for i = 1 : dims_with_vectors_and_scalars
      idx(i) = prepareIdx (varargin{i}, xsz(i));
    endfor

    ## If there are more varargin inputs than x dimensions, then input tests
    ## have verified that they are just scalars, so add [1 1 1 1 1 ... 1] to
    ## those dims to perform concatenation along those dims.
    if (n_rpts > n_xdims)
      for i = n_xdims + (1 : (n_rpts - n_xdims))
        idx(i) = ones (1, varargin{i});
      endfor
    endif

    ## Use completed idx to specify repetition of x values in all dimensions.
    ## The trailing ":" will take care of cases where n_xdims > n_rpts.
    retval = x(idx{:}, :);

  endif

endfunction



function x = repmatnew(A, m, varargin)

  if (nargin < 2)
    print_usage ();
  endif

  if (nargin == 3)
    n = varargin{1};
    if (! isempty (m) && isempty (n))
      m = m(:).';
      n = 1;
    elseif (isempty (m) && ! isempty (n))
      m = n(:).';
      n = 1;
    elseif (isempty (m) && isempty (n))
      m = n = 1;
    else
      if (all (size (m) > 1))
        m = m(:,1);
        if (numel (m) < 3)
          n = n(end);
        else
          n = [];
        endif
      endif
      if (all (size (n) > 1))
        n = n(:,1);
      endif
      m = m(:).';
      n = n(:).';
    endif
  else
    if (nargin > 3)
      ## input check for m and varargin
      if (isscalar (m) && all (cellfun ("numel", varargin) == 1))
        m = [m varargin{:}];
        n = [];
      else
        error ("repmat: all input arguments must be scalar");
      endif
    elseif (isempty (m))
      m = n = 1;
    elseif (isscalar (m))
      n = m;
    elseif (ndims (m) > 2)
      error ("repmat: M has more than 2 dimensions");
    elseif (all (size (m) > 1))
      m = m(:,1).';
      n = [];
    else
      m = m(:).';
      n = [];
    endif
  endif
  idx = [m, n];

  if (all (idx < 0))
    error ("repmat: invalid dimensions");
  else
    idx = max (idx, 0);
  endif

  if (numel (A) == 1)
    ## optimize the scalar fill case.
    if (any (idx == 0))
      x = resize (A, idx);
    else
      x(1:prod (idx)) = A;
      x = reshape (x, idx);
    endif
  elseif (ndims (A) == 2 && length (idx) < 3)
################################################################################
## See bug 59840 for reasons for different methods
    m=rows(A);
    n=columns(A);
    if (issparse(A)||((m+n+sum(idx)<175)&&(isnumeric(A)||islogical(A))))
      x=kron(ones(idx),A);
    else
      ## indexing is now faster, so we use it rather than kron.
      p = idx(1); q = idx(2);
      x = reshape (A, m, 1, n, 1);
      x = x(:, ones (1, p), :, ones (1, q));
      x = reshape (x, m*p, n*q);
    endif
################################################################################
  else
    aidx = size (A);
    ## ensure matching size
    idx(end+1:length (aidx)) = 1;
    aidx(end+1:length (idx)) = 1;
    ## create subscript array
    cidx = cell (2, length (aidx));
    for i = 1:length (aidx)
      cidx{1,i} = ':';
      cidx{2,i} = ones (1, idx (i));
    endfor
    aaidx = aidx;
    ## add singleton dims
    aaidx(2,:) = 1;
    A = reshape (A, aaidx(:));
    x = reshape (A (cidx{:}), idx .* aidx);
  endif

endfunction

function retval = repelemnew(x, varargin)

  if (nargin <= 1)
    print_usage ();

  elseif (nargin == 2)

    R = varargin{1};

    if (isscalar (R))

      if (! isvector (x))
        error (["repelem: %dD Array requires %d or more input " ...
                "arguments, but only %d given"], ...
               ndims (x), ndims (x) + 1, nargin);
      endif

      if (iscolumn (x))
        ## element values repeated R times in a col vector
        retval = x.'(ones (R, 1), :)(:);
      else
        ## element values repeated R times in a row vector
        retval = x(ones (R, 1), :)(:).';
      endif

    elseif (isvector (x) && isvector (R))

      ## vector x with vector repeat.
      if (numel (R) != numel (x))
        error (["repelem: R1 must either be scalar or have the same " ...
                "number of elements as the vector to be replicated"]);
      endif

      ## Basic run-length decoding in function prepareIdx returns
      ## idx2 as a row vector of element indices in the right positions.
      idx2 = prepareIdx (R);
      ## Fill with element values, direction matches element.
      retval = x(idx2);

    else # catch any arrays passed to x or varargin with nargin==2
      error (["repelem: when called with only two inputs they must be " ...
              "either scalars or vectors, not %s and %s."],
             typeinfo (x), typeinfo (R));
    endif

  elseif (nargin == 3)  # special optimized case for 2-D (matrices)

    ## Input Validation
    xsz = size (x);
    vector_r = ! (cellfun (@numel, varargin) == 1);

    ## 1. Check that all varargin are either scalars or vectors, not arrays.
    ##    isvector gives true for scalars.
    if (! (isvector (varargin{1}) && (isvector (varargin{2}))))
      error ("repelem: R1 and R2 must be scalars or vectors");

    ## 2. check that any repeat vectors have the right length.
    elseif (any (cellfun (@numel, varargin(vector_r)) != xsz(vector_r)))
      error (["repelem: R_j vectors must have the same number of elements " ...
              "as the size of dimension j of X"]);
    endif

    ## Create index arrays to pass to element.
    ## (It is no slower passing to prepareIdx
    ## than checking and doing scalars directly.)
################################################################################
## See bug 59840 for reasons for different methods
    if (isnumeric(x)||islogical(x))&&ismatrix(x)&&(issparse(x)||sum([xsz varargin{1} varargin{2}])<175)
      retval=kron(x,ones(varargin{1},varargin{2}));
    else
      idx1 = prepareIdx (varargin{1}, xsz(1));
      idx2 = prepareIdx (varargin{2}, xsz(2));

    ## The ":" at the end takes care of any x dimensions > 2.
      if issparse()
        retval = x(idx1, idx2);
      else
        retval = x(idx1, idx2, :);
      endif
    end
##################################################################################
  else  # (nargin > 3)

    ## Input Validation
    xsz = size (x);
    n_xdims = numel (xsz);
    vector_r = ! (cellfun (@numel, varargin) == 1);

    ## 1. Check that all repeats are scalars or vectors
    ##    (isvector gives true for scalars);
    if (! all (cellfun (@isvector, varargin(vector_r))))
      error ("repelem: R_j must all be scalars or vectors");

    ## 2. Catch any vectors thrown at trailing singletons,
    ##    which should only have scalars;
    elseif (find (vector_r, 1, "last") > n_xdims)
      error ("repelem: R_j for trailing singleton dimensions must be scalar");

    ## 3. Check that the ones that are vectors have the right length.
    elseif (any (cellfun (@numel, varargin(vector_r)) != xsz(vector_r)))
      error (["repelem: R_j vectors must have the same number of elements " ...
              "as the size of dimension j of X"]);

    endif

    n_rpts = nargin - 1;
    dims_with_vectors_and_scalars = min (n_xdims, n_rpts);

    ## Preallocate idx which will contain index array to be put into element.
    idx = cell (1, n_rpts);

    ## Use prepareIdx() to fill indices for dimensions that could be
    ## a scalar or a vector.
    for i = 1 : dims_with_vectors_and_scalars
      idx(i) = prepareIdx (varargin{i}, xsz(i));
    endfor

    ## If there are more varargin inputs than x dimensions, then input tests
    ## have verified that they are just scalars, so add [1 1 1 1 1 ... 1] to
    ## those dims to perform concatenation along those dims.
    if (n_rpts > n_xdims)
      for i = n_xdims + (1 : (n_rpts - n_xdims))
        idx(i) = ones (1, varargin{i});
      endfor
    endif

    ## Use completed idx to specify repetition of x values in all dimensions.
    ## The trailing ":" will take care of cases where n_xdims > n_rpts.
    retval = x(idx{:}, :);

  endif

endfunction

## Return a row vector of indices prepared for replicating.
function idx = prepareIdx (v, n)

  if (isscalar (v))
    ## will always return row vector
    idx = [1:n](ones (v, 1), :)(:).';

  else
    ## This works for a row or column vector.

    ## Get ending position for each element item.
    idx_temp = cumsum (v);

    ## Set starting position of each element to 1.
    idx(idx_temp + 1) = 1;

    ## Set starting position of each element to 1.
    idx(1) = 1;

    ## Row vector with proper length for output
    idx = idx(1:idx_temp(end));

    ## with prepared index
    idx = (find (v != 0))(cumsum (idx));

  endif

endfunction


(file #50689, file #50690, file #50691, file #50692)

Anonymous
Fri 08 Jan 2021 06:37:19 PM UTC, comment #3: 

Last word on the for loop, it doesn't calculate the same value for x as the other solutions.  It doesn't set the last row to anything so it is all zeros.  Since it isn't correct, it's unfair to compare on performance.

Rik <rik5>
Group administrator
Fri 08 Jan 2021 06:14:40 PM UTC, comment #2: 

A true replacement for repmat needs to not only accept an arbitrary number of dimensions as pointed out in comment #1, but it also needs to accept mostly arbitrary input including non-numeric input.

This code, for example, is perfectly acceptable


x.a = "123";
x.b = 456;
y = repmat (x, [2, 3])


but fails with the proposed repmatrixfast.

Also, there needs to be input validation to check the elements are correct.  Your implementation is the core of the algorithm, but I suspect that a large fraction of the difference is simply that required validation.  This could be checked with tic/toc calls around just the input validation section.

The beauty of Free Software is that you can go look at the source  code yourself.  The repmat.m function already has a specialization for 2-D inputs.  The code is


  elseif (ndims (A) == 2 && length (idx) < 3)
    if (issparse (A))
      x = kron (ones (idx), A);
    else
      ## indexing is now faster, so we use it rather than kron.
      m = rows (A); n = columns (A);
      p = idx(1); q = idx(2);
      x = reshape (A, m, 1, n, 1);
      x = x(:, ones (1, p), :, ones (1, q));
      x = reshape (x, m*p, n*q);
    endif


Form the comments, some programmer has checked the relative speed of indexing versus kron and decided on indexing.  If I excerpt just this piece of code in to a new function


function x = repmatidx (A, p, q)
  ## indexing is now faster, so we use it rather than kron.
  m = rows (A); n = columns (A);
  #p = idx(1); q = idx(2);
  x = reshape (A, m, 1, n, 1);
  x = x(:, ones (1, p), :, ones (1, q));
  x = reshape (x, m*p, n*q);
endfunction


and then add it in to the test script I get these timings


Origional attempt
Elapsed time is 2.06363 seconds.
with faster repmat
Elapsed time is 0.880177 seconds.
with faster repmatidx
Elapsed time is 1.33363 seconds.
single repmat call
Elapsed time is 1.28122 seconds.
with faster repmat
Elapsed time is 0.560874 seconds.
loop solution
Elapsed time is 0.278732 seconds.


So the kron solution is still faster.  Note that in repmat.m there is a special case path for sparse matrices.  One could imagine changing the test from "issparse (x)" to "isnumeric (x)" to get all numeric matrices.  I tried a new function


function x = repmatkron (A, p, q)

  x = kron (ones ([p, q]), A);

endfunction


which will only work for numeric matrices.  When I run the original test code the runtime is 0.70 versus the 0.88 of your implementation.  The difference is probably due to the minor amount of input argument decoding your function has.  Since kron() works on logical matrices (but not character matrices) it might be enough to change the special case to


if (isnumeric (A) || islogical (A))


Separately, the reason why the 'loop' scenario is fast is because it isn't really a loop.  The top of the script defines 'i = 2' and thus the loop


for k=[1:i-1]
===
for k=[1:2-1]
===
for k=[1:1]


So that loop just executes once.  In the second loop, nv =4, and so it isn't much of a loop either


for k=i+1:nv
===
for k=2+1:4
===
for k=3:4


So this loop executes twice for the values '3' and '4'.
-verbatim-

Rik <rik5>
Group administrator
Fri 08 Jan 2021 01:16:50 PM UTC, comment #1: 

Your timing data is acknowledged. It appears from your code that you are only looking for 2-dimensional replication, not for 3 or more dimensions. Is this correct? If that is the case, then that could account for speed differences. The standard repmat and repelem use varargin to allow for multiple arguments to maintain Matlab compatibility for calls like repmat(rand(2,3),2,3,4,5).

Would you be willing to write an improved repmat and repelem that allows for multiple arguments while maintaining the speed of your kron-based approach?

Anonymous
Fri 08 Jan 2021 01:49:29 AM UTC, original submission:  

repmat and repelem are slower than they need to be and in one case demonstrated they are slower than using loops.

Both repmat and repelem are slower than using kron with numeric matrix input and repeating in only one direction, repmat is also slower when repeating in both directions. Could my kron implementation be incorporated into repmat and repelem where it is faster? Also, why is my loop implementation fastest in the last example?


function repmatrepelemfaster
q=reshape(1:6,[],2)

repmat(q,2)
repmatmatrixfast(q,2)

repelem(q,2,3)
repelemmatrixfast(q,2,3)

ab=[1 500;
500 1;
50 50]

for k=1:size(ab,1)
  a=ab(k,1);
  b=ab(k,2);
  fprintf('b=%d, c=%d\n',a,b)
  disp('inbuilt repmat');
  tic
  for k=1:1000
    repmat(q,a,b);
  end
  toc
  disp('repmat using kron')
  tic
  for k=1:1000
    repmatmatrixfast(q,a,b);
  end
  toc
  disp('inbuilt repelem')
  tic
  for k=1:1000
    repelem(q,a,b);
  end
  toc
  disp('repelem using kron')
  tic
  for k=1:1000
    repelemmatrixfast(q,a,b);
  end
  toc
end

fprintf('What I was doing when I observed the issue:\nInsert a row into a column vector repeting the elements of the vector\nThis was used to plot a slice through a multivariable function.\nSpeed became an issue when I tried to intergrate the function.\n');

np=500;
xi=rand(1,np);
xo=rand(4,1);
i=2;
disp('Origional attempt')
tic
for l=1:10000
  x=[repmat(xo(1:i-1),1,np);xi;repmat(xo(i:end),1,np)];
end
toc
disp('with faster repmat')
tic
for l=1:10000
  x=[repmatmatrixfast(xo(1:i-1),1,np);xi;repmatmatrixfast(xo(i:end),1,np)];
end
toc
disp('single repmat call')
tic
for l=1:10000
  x=zeros(size(xo,1)+1,size(xi,2));
  x([1:i-1 i+1:end],:)=repmat(xo,1,np);
  x(i,:)=xi;
end
toc
disp('with faster repmat')
tic
for l=1:10000
  x=zeros(size(xo,1)+1,size(xi,2));
  x([1:i-1 i+1:end],:)=repmatmatrixfast(xo,1,np);
  x(i,:)=xi;
end
toc
disp('loop solution')
tic
for l=1:10000
  nv=size(xo,1);
  x=zeros(size(xo,1)+1,size(xi,2));
  x(i,:)=xi;
  for k=[1:i-1]
    x(k,:)=xo(k);
  end
  for k=i+1:nv
    x(k,:)=xo(k-1);
  end
end
toc

end
function r=repmatmatrixfast(a,b,c)
if numel(b)==2
  c=b(2);
  b=b(1);
elseif nargin()==2
  c=b;
end
r=kron(ones(b,1),a);
r=kron(ones(1,c),r);
end
function r=repelemmatrixfast(a,b,c)
r=kron(a,ones(b,1));
r=kron(r,ones(1,c));
end


Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #50748:  repmatrepelemfaster.m added by None (21KiB - text/plain)
file #50743:  repmattimekronindx added by None (2MiB - application/octet-stream)
file #50744:  repmatrepelemfaster.m added by None (21KiB - text/plain)
file #50727:  objf.m added by rik5 (321B - text/x-matlab)
file #50692:  repmatoldnewsametime2.eps added by None (221KiB - application/postscript)
file #50689:  repelemoldnewsametime.eps added by None (220KiB - application/postscript)
file #50690:  repelemoldnewsametime2.eps added by None (220KiB - application/postscript)
file #50691:  repmatoldnewsametime.eps added by None (222KiB - application/postscript)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by nrjank
  • -email is unavailable- added by rik5 (Posted a comment)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 10 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-09-23 nrjank Carbon-Copy- Added -email is unavailable-
    2021-01-19 None Attached File- Added repmatrepelemfaster.m, #50748
    2021-01-18 None Attached File- Added repmattimekronindx, #50743
        Attached File- Added repmatrepelemfaster.m, #50744
    2021-01-17 rik5 Attached File- Added objf.m, #50727
        StatusNone In Progress
    2021-01-10 None Attached File- Added repmatoldnewsametime.eps, #50691
        Attached File- Added repmatoldnewsametime2.eps, #50692
    2021-01-10 None Attached File- Added repelemoldnewsametime.eps, #50689
        Attached File- Added repelemoldnewsametime2.eps, #50690

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code