How to implement Outer with loops

Question

I would like to implement a simplified version of Outer (Outer[f, vec1, vec2,..]) through a looping construct.

Here are is one thing that works, but is not quite down to the level I'd like it to be on (think, say, what is a available in C). Using Table one could do the following:

v1 = {a, b, c};
v2 = {A, B, C, D};
Table[Table[f[i, j], {j, v2}], {i, v1}] === Outer[f, v1, v2]
True

I know that the result should be of the form (this is correct for n vectors):

lsts = {v1, v2};
res = ConstantArray[0, Join[Length /@ lsts, {Length[lsts]}]]

A Do loop would look like this:

Do[
 Do[
  res[[j, i, 1]] = lsts[[1, j]];
  res[[j, i, 2]] = lsts[[2, i]];
  , {i, Length[lsts[[2]]]}]
 , {j, Length[lsts[[1]]]}]

res == Outer[List, lsts[[1]], lsts[[2]]]
True

My question is how would I generalize the do loop for n vectors?

Mathics does this via recursion. Looking at the output of a compiled outer did not give an idea either:

cf = Compile[{}, Outer[List, {1, 2, 3}, {4, 5, 6}]];
Needs["CompiledFunctionTools`"]
CompilePrint[cf]

Why would you like to do this? If this is a quest for better understanding you might want to skip loops altogether and if you are trying to improve performance you ought to look elsewhere since loops are not efficient in Mathematica. — Sascha, Dec 10 '16 at 23:26
You can build your Do-loop code using metaprogramming. However one usually wishes to go the other direction in Mathematica, i.e. convert verbose loops into concise high-level abstractions such as Outer, therefore I echo Sascha's question: why? If you explain the purpose behind this question it is likely that you will get more useful replies that if you do not. — Mr.Wizard, Dec 11 '16 at 05:41
@Sascha, I'd like to implement a Outer in C for integer packed arrays. — TheGuest, Dec 11 '16 at 09:07
@Mr.Wizard, while it's true that in most cases vectorization is the way to go in mathematica it's by no means always the case. For example the use of loops in Compile can outperform vectorized operations, esp. when the compilation target is C. Maybe this forum is not quite the right place to ask, but in a C forum people might have a hard time to understand Outer. — TheGuest, Dec 11 '16 at 09:09
If you're going for performance in a particular situation, I would suggest amending your question with a simple example problem for benchmarking, and maybe an algorithm to beat. People go ape for questions like that around here... — Simon Rochester, Dec 11 '16 at 10:47
I'm not clear on the goal here. Are you using Mathematica to prototype an algorithm which you will write in C, or are you after a CompiledFunction in Mathematica which you can use in place of Outer? — Simon Woods, Dec 11 '16 at 11:58
My answer to "How to implement the general array broadcasting method from NumPy?" question contains "C" implementation of broadcasting, which "overlaps" with Outer. — jkuczm, Jan 21 '17 at 14:50

score 1 · Answer 1 · answered Jan 19 '17 at 16:20

If you want an implementation in C, this might be quite close. Here's an example with three vectors:

{x1, x2, x3} = {{1, 2, 3}, {1, 2, 3, 4}, {1, 2}};

out = ConstantArray[0, {Length[x1], Length[x2], Length[x3]}];
For[i1 = 1, i1 <= Length[x1], i1++,
 For[i2 = 1, i2 <= Length[x2], i2++,
  For[i3 = 1, i3 <= Length[x3], i3++,
   out[[i1, i2, i3]] = f[x1[[i1]], x2[[i2]], x3[[i3]]]
   ]
  ]
 ]

out == Outer[f, Range[3], Range[4], Range[2]]
(* True *)

Here's a List example:

{x[1], x[2], x[3]} = {{1, 2, 3}, {1, 2, 3, 4}, {1, 2}};
nArgs = 3;
out = ConstantArray[
   0, {Length[x[1]], Length[x[2]], Length[x[3]], nArgs}];
For[i[1] = 1, i[1] <= Length[x[1]], i[1]++,
 For[i[2] = 1, i[2] <= Length[x[2]], i[2]++,
  For[i[3] = 1, i[3] <= Length[x[3]], i[3]++,
   For[i[0] = 1, i[0] <= nArgs, i[0]++,
    out[[i[1], i[2], i[3], i[0]]] = x[ i[0] ][[ i[ i[0] ] ]]
    ]
   ]
  ]
 ]

out == Outer[List, Range[3], Range[4], Range[2]]
(* True *)

Notice, how I switched over to indexed variables x[1], x[2], etc. to make generalization easier. i[0] here refers to the position of the argument.

This does indeed seem like a task for recursion, but I don't have time to go that far.

How to implement Outer with loops

1 Answers1