Asked by Adi gahlawat
on 27 Jul 2013

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

*No products are associated with this question.*

Answer by Azzi Abdelmalek
on 27 Jul 2013

Edited by Azzi Abdelmalek
on 27 Jul 2013

Accepted answer

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

Show 4 older comments

Azzi Abdelmalek
on 27 Jul 2013

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

Answer by Roger Stafford
on 27 Jul 2013

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
on 27 Jul 2013

Hi Roger,

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

Adi

Answer by Azzi Abdelmalek
on 28 Jul 2013

Edited by Azzi Abdelmalek
on 28 Jul 2013

**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

Answer by Richard Brown
on 29 Jul 2013

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

Show 4 older comments

Azzi Abdelmalek
on 29 Jul 2013

Almost, the same result

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

Richard Brown
on 29 Jul 2013

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 by Jan Simon
on 29 Jul 2013

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
on 29 Jul 2013

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

Answer by Richard Brown
on 29 Jul 2013

Edited by Richard Brown
on 29 Jul 2013

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
on 29 Jul 2013

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

Richard Brown
on 29 Jul 2013

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 by Andrei Bobrov
on 27 Jul 2013

Edited by Andrei Bobrov
on 30 Jul 2013

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));

Opportunities for recent engineering grads.

## 0 Comments