How to avoid for loops when generating index arrays?

I often find myself coding nested for loops to generate vectors of integer indices. For example:

n = 4;
i = 1;
for L = 0:n
    for M = -L:L
        l(i) = L;
        m(i) = M;
        i = i+1;
    end
end

All I need are the vectors "l" and "m". I can preallocate to save some speed, but my real problem is having to use the for loops as sometimes the index vectors I need to create have many more nested for loops whose (note that the inner loop index depends on the outer loop index).

Is there a simple way to avoid using loops to generate index vectors like these?

Answer 1

[Matt J]

Here's another method, less memory consuming than NDGRID

    mm=sparse(-n:n);
    ll=sparse(0:n);

    map=bsxfun(@le,abs(mm.'),ll);

    idx=nonzeros(bsxfun(@times,map,1:length(ll)  ));
       l=full(ll(idx));

    idx=nonzeros(bsxfun(@times,map,(1:length(mm)).'))  ;
       m=full(mm(idx));

Answer 2

[Roger Stafford]

For your particular problem you can do this:

 M = (0:n*(n+2))';
 L = floor(sqrt(M));
 M = M-L.*(L+1);

(I've used uppercase letters, 'L' and 'M', in place of your lowercase 'l' and 'm'.)

As with Matt Kindig, I am not sure this will be any faster than your for-loops. Time it with a large value for n and see.

Answer 3

[Matt J]

Here's a way to do it using NDGRID. It's not apriori obvious whether for loops would or would not be faster. It depends what you plan to reuse.

[mg,lg]=ndgrid(-n:n,0:n);

idx=abs(mg)<=lg;
 l=lg(idx).',
 m=mg(idx).',

[Oliver]

Thanks Matt, this is definitely a simple way to do it. However, it requires creation of an intermediate array which can be much larger than the final array. Is there a solution which is efficient in terms of both memory and vectorization?

[Matt J]

To answer that meaningfully, you need to give us more of a sense of the scope of what you're doing. I cannot see why the intermediate array "can be much larger than the final array". In this example, they are less than twice its size. As I also mentioned, it's not clear what you might be able to reuse. If you plan to generate further arrays, it might be worth keeping lg and mg around, regardless of the memory they consume.

[Oliver]

Consider the case of n = 1. In this case the vectors l and m should be:

l = [0, 1, 1, 1]; m = [0, -1, 0, 1];

Using ndgrid would create lg and mg which are each 2-by-3 arrays, clearly larger than the final result in which l and m are 1-by-4 arrays. If you consider something much larger, such as n = 10 you get size(lg) = [21,11], whereas size(l) = [1 121]. In this case the intermediate array lg uses nearly twice as much memory as the final array l.

[Matt J]

In this case the intermediate array lg uses nearly twice as much memory as the final array l.

Again, why do we care? Even for n=300, lg and mg consume a measly 2MB. And as I've said (3 times now), you might be able to reuse lg and mg to generate other stuff.

[Oliver]

The reason that I care is that I have only given you a simplified example. In practice I have things like:

for N = 0:2:30
    for L = 0:N
        for M = -L:L
            for P = 0:M
                for Q = -P:P
                    ...
                end
            end
        end
    end
end

So making the intermediate arrays in this case would require a 5-dimensional array that takes up a lot of space.

[Matt J]

I'm starting to think Sean's advice about sticking with for-loops is the best one. There can definitely be ways to cut down on the loop nesting (see my newest Answer based on cell arrays), but the required form would depend on the body of the original set of for-loops.

Answer 4

[Sean de Wolski]

doc meshgrid
doc ndgrid %?

:)

And of course, depending on your application, two nested for-loops or bsxfun() might be better.

[Oliver]

I can't use ndgrid or meshgrid because the inner loop depends on the outer loop. If i had something like

i = 1;
for L = 3:6
    for M = 1:5
        l(i) = L;
        m(i) = M;
        i = i+1;
    end
end

then I could use:

[l,m] = meshgrid(3:6,1:5);

But, this won't work with the dependency.

[Sean de Wolski]

Just use the for-loops, they'll be the fastest by far. If you want to disguise it, write a function that takes L and M and returns l and m.

cellfun and arrayfun are slow and converting between cells and numeric types is slow. The above with preallocation will be pretty quick.

Answer 5

[Matt Kindig]

It's kind of hack-y, but it gives the same output as your original posting:

    n=4;
    l = cell2mat(arrayfun(@(x) x*ones(1,2*x+1), 0:n, 'uni', false));
    m= cell2mat( arrayfun(@(x) (-x:1:x), 0:n, 'uni', false));

Keep in mind that this may very well be slower than for-loops--I haven't done any timing comparisons.

Answer 6

[Matt J]

Here's a way to eliminate one nested loop

    l=cell(1,n+1);
    m=l;
    for L=0:n

        i=L+1;
        m{i}=-L:L;
        l{i}=m{i};
         l{i}(:)=L;

    end

    l=[l{:}],
    m=[m{:}],

[Sean de Wolski]

I'd be surprised if this is faster due to the cell array conversions. I guess one of us will have to run a timing test.

[Matt J]

For n=1000 I get this,

Original Approach:

Elapsed time is 0.093130 seconds.

Cell-Based Approach

Elapsed time is 0.027393 seconds.

I think the vectorization inherent in

        m{i}=-L:L;
        l{i}=m{i};
         l{i}(:)=L;

trumps the overhead from the cell conversion.