Generating a particular sequnce of numbers

Hi,

given a variable natural number d, I'm trying to generate a sequence of the form:

[1 2 1 3 2 1 4 3 2 1.......d d-1 d-2......3 2 1].

I don't want to use for loop for this process, does anyone know a better (faster) method. I tried the colon operator without any success.

Thank you.

Adi

Answer 1

[Azzi Abdelmalek]

d=4
cell2mat(arrayfun(@(x) x:-1:1,1:d,'un',0))

[Adi gahlawat]

Thanks a lot Azzi, I'll give it a shot. At least now I have a method to compare with the for loop.

[Youssef Khmou]

that is fast method, i voted for it, the only thing that is left is explaining why the option 'UniformOutput' is set to false, (just for further questions including arrayfun)

[Adi gahlawat]

Hi Youssef, thanks or the input. The method proposed by Azzi is much slower than utilising for loops. Check out the following for d=df+1:

df=3;
tic
A1=cell2mat(arrayfun(@(x) x:-1:1,1:df+1,'un',0))';
toc
tic
A2=zeros(df+1,df+1);
for i=1:df+1
  A2(i,i:df+1)=1:df+2-i;
end
A2=A2(:); % Converting matrix to a vector
A2=A2(A2~=0); % Removing zeros
toc

The output is:

Elapsed time is 0.007756 seconds.
Elapsed time is 0.000032 seconds.

As you can see, the for loop based method is order 2 faster.

[Azzi Abdelmalek]

Youcef, the result for each cell is a different size, that's why we set "Uniform Output" to False

[Azzi Abdelmalek]

Adi gahlawat, you can not compare by puting df=3, look at this:

df=1000;
tic
for k=1:500
A1=cell2mat(arrayfun(@(x) x:-1:1,1:df+1,'un',0))';
end
toc
tic
for k=1:500
A2=zeros(df+1,df+1);
for i=1:df+1
A2(i,i:df+1)=1:df+2-i;
end
A2=A2(:); % Converting matrix to a vector
A2=A2(A2~=0); % Removing zeros
end
toc
tic
for k=1:500
   N = df*(df+1)/2;
A = zeros(1,N);
n = 1:df;
A((n.^2-n+2)/2) = n;
A = cumsum(A)-(1:N)+1;
end
toc

Elapsed time is 5.617201 seconds.  % Azzi's answer
Elapsed time is 10.643052 seconds. % Adi's answer
Elapsed time is 4.755004 seconds.  % Stafford's answer

[Youssef Khmou]

ok thank you for the explanation .

[Jan]

Azzi's for loop approach is faster.

Answer 2

[Roger Stafford]

Here's another method to try:

 N = d*(d+1)/2;
 A = zeros(1,N);
 n = 1:d;
 A((n.^2-n+2)/2) = n;
 A = cumsum(A)-(1:N)+1;

[Adi gahlawat]

Hi Roger,

your method is excellent. It's about 2 times faster than my for loop based code. Much obliged.

Adi

Answer 3

[Azzi Abdelmalek]

Edit

This is twice faster then Stafford's answer

A4=zeros(1,d*(d+1)/2); % Pre-allocate
c=0;
for k=1:d
  A4(c+1:c+k)=k:-1:1;
  c=c+k;
end

[Jan]

Yes, this is exactly the kind of simplicity, which runs fast. While the one-liners with anonymous functions processed by cellfun or arrayfun look sophisticated, such basic loops hit the point. +1

I'd replace sum(1:d) by: d*(d+1)/2 . Anbd you can omit idx.

Answer 4

[Richard Brown]

Even faster:

k = 1;
n = d*(d+1)/2;
out = zeros(n, 1);

for i = 1:d
    for j = i:-1:1
        out(k) = j;
        k = k + 1;
    end
end

[Jan]

In my measurements this is remarkably slower than Azzi's loop approach. Which Matlab version are you using?

[Richard Brown]

2012a on 64 bit Linux. Can try R2013a tomorrow for interest's sake :-)

[Azzi Abdelmalek]

Test with d=3000 and for loop (500)

d=3000;
tic
for k=1:500
k = 1;
n = d*(d+1)/2;
out = zeros(n, 1);
for i = 1:d
  for j = i:-1:1
      out(k) = j;
      k = k + 1;
  end
end
end
toc
tic
for k=1:500
  A4=zeros(1,d*(d+1)/2); % Pre-allocate
c=0;
for k=1:d
A4(c+1:c+k)=k:-1:1;
c=c+k;
end
end
toc
isequal(A4,out')

Elapsed time is 23.695618 seconds. Richard's answer

Elapsed time is 17.408498 seconds. Azzi's answer

[Richard Brown]

Can you try again with a different loop counter name? I'm not in front of a computer right now so can't check whether that interferes or not...

[Azzi Abdelmalek]

Almost, the same result

Elapsed time is 22.940850 seconds.
Elapsed time is 16.967270 seconds.

[Jan]

Under R2011b I get for d=1000 and 500 repetitions:

 Elapsed time is 3.466296 seconds.  Azzi's loop
 Elapsed time is 3.765340 seconds.  Richard's double loop
 Elapsed time is 1.897343 seconds.  C-Mex (see my answer)

[Richard Brown]

I checked again, and I agree with Azzi. My method was running faster because of another case I had in between his and mine. The JIT was doing some kind of unanticipated optimisation between cases.

I get similar orders of magnitude results to Azzi for R2012a if I remove that case, and if I run in R2013a (Linux), his method is twice as fast.

Shame, I like it when JIT brings performance of completely naive loops up to vectorised speed :)

Answer 5

[Jan]

An finally the C-Mex:

#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[]) {
  mwSize d, i, j;
  double *r;

    d       = (mwSize) mxGetScalar(prhs[0]);
    plhs[0] = mxCreateDoubleMatrix(1, d * (d + 1) / 2, mxREAL);
    r       = mxGetPr(plhs[0]);
    for (i = 1; i <= d; i++) {
       for (j = i; j != 0; *r++ = j--) ;
    }
  }

And if your number d can be limited to 65535, the times shrink from 1.9 to 0.34 seconds:

#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[]) {
  uint16_T d, i, j, *r;

    d       = (uint16_T) mxGetScalar(prhs[0]);
    plhs[0] = mxCreateNumericMatrix(1, d * (d + 1) / 2, mxUINT16_CLASS, mxREAL);
    r       = (uint16_T *) mxGetData(plhs[0]);
    for (i = 1; i <= d; i++) {
       for (j = i; j != 0; *r++ = j--) ;
    }
  }

For UINT32 0.89 seconds are required.

[Richard Brown]

Nice. I imagine d would be limited to less than 65535, that's a pretty huge vector otherwise

Answer 6

[Richard Brown]

Also comparable, but not (quite) faster

n = 1:(d*(d+1)/2);
a = ceil(0.5*(-1 + sqrt(1 + 8*n)));
out = a.*(a + 1)/2 - n + 1;

[Richard Brown]

Potentially suffers from floating point errors, but I checked it up to d = 10000 :)

[Jan]

@Richard: How did you find this formula?

[Richard Brown]

If you look at the sequence, and add 0, 1, 2, 3, 4 ... you get

 n:  1  2  3  4  5  6  7  8  9 10
     1  3  3  6  6  6 10 10 10 10

Note that these are the triangular numbers, and that the triangular numbers 1, 3, 6, 10 appear in their corresponding positions, The a-th triangular number is given by

n = a (a + 1) / 2

So if you solve this quadratic for a where n is a triangular number, you get the index of the triangular number. If you do this for a value of n in between two triangular numbers, you can round this up, and invert the formula to get the nearest triangular number above (which is what the sequence is). Finally, you just subtract the sequence 0, 1, 2, ... to recover the original one.

Answer 7

[Andrei Bobrov]

out = nonzeros(triu(toeplitz(1:d)));

or

out = bsxfun(@minus,1:d,(0:d-1)');
out = out(out>0);

or

z = 1:d;
z2 = cumsum(z);
z1 = z2 - z + 1;
for jj = d:-1:1
    out(z1(jj):z2(jj)) = jj:-1:1;
end

or

out = ones(d*(d+1)/2,1);
ii = cumsum(d:-1:1) - (d:-1:1) + 1;
out(ii(2:end)) = 1-d : -1;
out = flipud(cumsum(out));