Alternatives to procedural loops and iterating over lists in Mathematica

Question

While there are some cases where a For loop might be reasonable, it's a general mantra – one I subscribe to myself – that "if you are using a For loop in Mathematica, you are probably doing it wrong". But For and Do loops are familiar to people who already know other programming languages. It is common for new Mathematica users to rely on loop constructs, only to become frustrated with performance or code complexity.

My question is: are there some general rules of thumb that would help new users select the appropriate alternative to a procedural For or Do loop, and are there other useful Mathematica constructs, aside from RandomVariate/Nest/Fold/Inner/Outer/Tuples that they should know about in order to avoid those loops?

---

In particular, there are a number of ways of iterating over lists:

1. The For loop:

   For[i = 1, i <= 10, i++, list[[i]] = func[ list[[i]]]
   

2. The Table function operating over each part of a list in turn:

   Table[func[ list[[i]] ], {i, Length[list]}]
   

3. The Do loop operating over each part of a list in turn:

   lst = {88, 42, 25, 75, 35, 97, 12};
   t = 9;
   Do[
     x = lst[[i]];
     t += Mod[x, t],
     {i, 1, Length[lst]}
   ];
   t
   

4. Mapping a function onto each element of a list:

   (3 - #)/(7 * #) & /@ list
   

5. And for nested lists, there are constructs like MapThread:

   f = Mod[#, #2] Floor[#3/#] &;
   
   a = {{18, 85, 22, 20, 39}, {17, 67, 76, 96, 58}, {40, 97, 56, 60, 53}};
   
   MapThread[f, a]
   

New users of Mathematica will usually pick options 1 or 3 and use C-style or Matlab approaches that involve the following multistep process:

1. Defining an empty list
2. Setting up a loop
3. Optionally define a variable (usually not localised) to equal the iterator, which is then used in subsequent calculations within the loop
4. Within each loop iteration, use the local variable defined to equal the iterator to redefine each element of that list in turn according to some function
5. If the list is multidimensional, nest another loop inside that.

What are some useful guides to help users coming from other languages to Mathematica to improve the conciseness and efficiency of their code by avoiding unnecessary loops?

Answer 1

Let's improve the specific cases given.

---

Case #1

Explicit loops are often counterproductive in Mathematica, not only taking more keystrokes, but also more execution time. They are also, in my opinion, more prone to mistakes.

Better ways are to use Do, Scan, or Map.

Do and Scan are (typically) appropriate for operations that do not accumulate a list of results, while Map and Table (a variant of Do) are (typically) used for ones that do. The third method uses Map (short form: /@) and accumulates a useless list of results; for this reason it is likely to be less memory efficient than the first two.

1. Do[Print[i], {i, 10}]
   

2. Scan[Print, Range@10]
   

3. Print /@ Range@10;
   

Average Timing, using PrimeQ in place of Print for 10^6 iterations:

---

Case #2

Since version 6 this form of Table is not necessary! Instead, use:

Table[func[i], {i, list}]

One might ask why use Table here at all. Indeed, this simple example can be written:

func /@ list

which is preferred. However, there are more complex cases where Table is far more elegant than the alternatives.

Average Timing for 10^6 real numbers using Sin for func:

---

Case #3

This could be improved in the manner of case #2, but there is a far better method available: Fold. Any time you want to iterate through a list, using the result of the previous "loop" along the way, look to Fold, efficient in both syntax and computation.

Fold[# + Mod[#2, #] &, 9, lst]

Average Timing for 10^6 integers (big savings!):

---

Case #4

There is nothing intrinsically wrong with this. However, all of the operations within the pure function have the Attribute Listable. Therefore, this function can directly accept a list without using Map, and this will be considerably more efficient.

(3 - #)/(7 * #) & @ list

Average Timing for a list of 10^6 real numbers:

---

Case #5

This is similar to case #4, but a little more complicated. Once again, each subfunction in f is Listable but this is not being leveraged. One can use Apply to pass the sub-lists in a as arguments to f:

f @@ a

Suppose that not all of the functions are Listable. I create a dummy addition function g that only accepts integers, not lists. I then include this in f, and try to Apply it again:

ClearAll[g, f]
g[n_Integer, m_Integer] := n + m
f = Mod[#, #2] Floor[#3/#] * g[#2, #] &;

f @@ a

But the result is incorrect. Once could go back to MapThread, but a better way, when possible, is to make g handle lists, which will usually be faster on large sets. Here are two ways of doing that. Give g the Listable attribute, and Mathematica will automatically thread over lists:

ClearAll[g, f]
SetAttributes[g, Listable]
g[n_Integer, m_Integer] := n + m
f = Mod[#, #2] Floor[#3/#] * g[#2, #] &;

f @@ a

Or, if the automatic threading via Listable breaks your function in some way, manually:

ClearAll[g, f]
g[n_Integer, m_Integer] := n + m
g[n_List, m_List] := MapThread[Plus, {n, m}]
f = Mod[#, #2] Floor[#3/#] * g[#2, #] &;

f @@ a

Average Timing for lists of 10^6 integers:

Answer 2

There are many alternative ways to approach various programming problems that do not use loops and are more efficient (and concise) in Mathematica. Most of them execute faster, but even where they do not, they are faster to type: development time matters, too!

Here are some rules of thumb for easier programming and iterating on lists.

1. Most arithmetic operations and many other functions are Listable, meaning that they operate element by element without having to set up a loop or operate explicitly on each element

Addition, multiplication and other standard operations work element by element when they are given conformable vectors, matrices, tensors or generalised lists:

list1 = {a, b, c};
list2 = {e, f, g};
list1 + list2
(* {a + e, b + f, c + g} *)

Listable functions include standard arithmetic and trigonometric functions, Log, Exp, Mod, Abs and a huge range of exponential, Bessel-related and other special functions. The following code returns the whole list of Listable functions:

With[{names = Names["System`*"]}, 
 Pick[names, MemberQ[Attributes[#], Listable] & /@ names, True] ]

(see also the answer by Jagra)

2. In Mathematica, you do not need to declare an empty list and then fill each element in with the desired value one by one

For example, typical procedural programming will create an empty vector, and replace each element in turn to create the desired data.

vec = Table[0, {100}]; For[i = 1, i <= 100, vec[[i]] = RandomReal[]+2; i++]; vec

In Mathematica, one can use RandomVariate and its cousins directly, or perhaps a Table command. This generalizes to other more complex definitions inside the Table function, and one can use the fact that many arithmetic operations are Listable to avoid looping. (Yes, I know this could be written as RandomReal[{2, 3}, 100], but I wanted a simple example of vectorized arithmetic operations.)

vec = RandomReal[{0, 1}, 100]+2.;

vec = Table[RandomReal[], {100}]+2.;

3. If the i-th element of your list depends on the previous elements, you can construct the list using Nest, Fold, FixedPoint and their variants NestList, FoldList and FixedPointList

Another common use of For loops is when output i depends on output i-1. In Mathematica, this is a canonical use of Nest and Fold, and their NestList and FoldList counterparts that also return the intermediate results. For example, here is the For loop way to create autoregressive noise:

vec = Table[0, {100}]; For[i = 2, i <= 100, 
  vec[[i]] = vec[[i - 1]] + RandomReal[]; i++];

Of course, you have to start at i=2 to avoid adding the zeroth part, which is the Head, List, to the vector.

Here is the FoldList way:

FoldList[#1 + #2 &, 0, RandomReal[{0, 1}, 99] ]

4. You can combine lists in very general ways using Inner and Outer.

People also often use nested For loops when they want to build up matrices depending on elements of two different vectors. Here is where the generalised structure of Outer is useful:

Outer[N@Kurtosis[#1[#2]] &, {StudentTDistribution, ExponentialDistribution, 
  ChiSquareDistribution}, Range[5, 15] ]

(* {{9., 6., 5., 4.5, 4.2, 4., 3.85714, 3.75, 3.66667, 3.6, 3.54545}, {9., 9., 
  9., 9., 9., 9., 9., 9., 9., 9., 9.}, {5.4, 5., 4.71429, 4.5, 4.33333, 4.2, 
  4.09091, 4., 3.92308, 3.85714, 3.8}} *)

Another way to avoid nested loops is with Tuples, as discussed here. This nested For loop:

l = 0;
For[i = 0, i < 10, i++,
 For[j = 0, j < 10, j++,
  For[k = 0, k < 10, k++,
   Xarray[l] = A[i, j, k];   
   Print[Xarray[l]];
   l++;
   ]
  ]
 ]
Xarray[5]

Can be rewritten faster as:

Xarray = A @@@ Tuples[Range[0, 9], 3];

Answer 3

How about operations like arithmetic on lists, matrixes, & tensors or Dot and DotProduct?

list1 = {a, b, c};
list2 = {d, e, f};
matrix1 = {{a, b, c}, {d, e, f}};
matrix2 = {{a, b, c}, {d, e, f}, {g, h, i}};

tensor1 = {{{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, {{1, 2, 3}, {4, 5, 
     6}, {7, 8, 9}}, {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}}, {{{a, b, 
     c}, {d, e, f}, {g, h, i}}, {{a, b, c}, {d, e, f}, {g, h, 
     i}}, {{a, b, c}, {d, e, f}, {g, h, i}}}}

list1 + n
list1 - n
list1 * n
list1 / n
list1 + list2
list1 - list2
list1 * list2
list1 / list2

matrix1 . list1
matrix1 . matrix2
tensor1 . list1
tensor1 . matrix2
tensor1 . matrix2 . list1

Just a sampling of things, some of which would require some extensive procedural solutions.

"For" my part (pun intended), I stopped thinking in loops and started thinking functionally. More so, I began to think more like a mathematician than a programer. The better I do this the more efficiently I program in Mathematica.

It would help if one could tie someone's "loop" or "procedural" arm to their body, like they do with boxers who need to switch their dominant hand.

Ultimately, I think one has to give themselves over to the functional language. Maybe one can translate one language to another but poetry in English doesn't have the same meaning translated into Italian nor Italian to English. When one can dream in the language they'll have it.

---

Some additional ideas...

The Vector analysis package

Needs["VectorAnalysis`"]

Includes:

• DotProduct[],
• CrossProduct[], and
• ScalarTripleProduct[]

Brief examples from the documentation:

DotProduct[]

Dot product of two Cartesian vectors:

In2:= a = {1, 2, 5};

In3:= b = {2, 3, -7};

In[4]:= DotProduct[a, b]

Out[4]= -27

CrossProduct[]

Find the cross product of a pair of vectors:

In2:= a = {1, 3, 5};

In3:= b = {-4, 7, 1};

In[4]:= CrossProduct[a, b]

Out[4]= {-32, -21, 19}

ScalarTripleProduct[]

Find the equation of the plane passing through the points with position vectors r1, r2, and r3:

In[7]:= r = {x, y, z};

In[8]:= ScalarTripleProduct[r - r1, r2 - r1, r3 - r1] == 0

Out[8]= -111 + 68 x - 29 y - 22 z == 0

---

It just doesn't stop...

Consider one example of functional decompositions, the CholeskyDecomposition[].

Maybe this doesn't count in this specific discussion because it demonstrates a specific type of solution available for use in Mathematica, which essentially hides how it gets the solution from the user. On the other side, I'd argue that the functional paradigm in Mathematica has made these kinds of solutions and functionality more readily available. The Wikipedia entry cited above shows how procedural approaches attack the problem. One could replace every procedural construct in that code with a functional one and produce your own functional Cholesky Decomposition pretty readily. But Mathematica has already done that for us.

The other decompositions:

LUDecomposition[]
QRDecomposition[]
SingularValueDecomposition[]
SchurDecomposition[]
HessenbergDecomposition[]
JordanDecomposition[]...
(* I've likely missed a few *)

probably do the same thing.