Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question

migrated from stackoverflow.com Jul 5 '12 at 13:12

This question came from our site for professional and enthusiast programmers.

    
How about operations like arithmetic on lists, matrixes, & tensors or Dot and DotProduct? "For" my part (pun intended), I stopped thinking in loops and started thinking functionally. It would help if one could tie someones "loop" 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. –  Jagra Jul 5 '12 at 12:53
    
There are cases where I find for loops more readable. For example R = S (*list of elemetns*); Do[ R[[i]] = f[S[[i]], G[[i]]]; , {i, Length@S}]; vs R = f[#[[1]], #[[2]]] & /@ Transpose[{S, G}] The second one is more concise, but G[[i]] conveys the intent better then #[[1]]. –  Ajasja Jul 5 '12 at 12:54
4  
@Ajasja, even more concise and readable is R=Thread@f[S,G] –  Simon Woods Jul 5 '12 at 13:28
3  
Just one other thing one may want to add as a caveat: the preferred choice of construct can change radically when the goal is to use it with Compile. –  Jens Jul 5 '12 at 15:12
1  
@Jens yes, that seems to be ignored most of the time. readability is also dependent on the person: I have some code to transform a density matrix that looks like this: MapThread[Join,{rtmp,Rest /@ Flatten[Conjugate@rtmp, {{2}, {1}}]}]. Some may read this immediately, but I need a commented-out iterative version to be able to re-read it later (this one's faster though). –  acl Jul 5 '12 at 15:34

3 Answers 3

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:

{For, Do, Scan, Map} = {1.5, 0.469, 0.375, 0.375}


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:

{Old, New} = {0.1438, 0.1124}


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!):

{Do loop, Fold} = {2.843, 0.0968}


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:

{Map, Listable} = {0.2218, 0.0512}


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:

{MapThread, Apply} = {2.75, 1.016}

share|improve this answer
18  
Wow! That was quick! The OP must be amazed! ... Wait .. –  belisarius Apr 30 '11 at 0:30
2  
@belisarius lol er um yeah. BTW, notice case #2. –  Mr.Wizard Apr 30 '11 at 0:33
    
@Mr. I use #2 all the time .... but the version in the question. Habits are difficult to change (said the monk) –  belisarius Apr 30 '11 at 5:34
    
@Mr.Wizard I think it'd be nice to give examples of the Do and Map functions you use in Case #3 and #4, respectively. Lets people test out the difference in speed themselves! –  Vincent Tjeng Mar 31 '13 at 3:46
    
@Vincent Those were accidentally left out of the question for a time, but I restored them shortly after your comment. (Thanks.) Now almost a year later do you have any other suggestions for improving this Q&A? –  Mr.Wizard Jan 4 at 21:43

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];
share|improve this answer
    
Great description. Thanks. And, in the last example, for loops nested, Xarray[l] would be Xarray[[l]], right? –  d555 Sep 8 '13 at 21:44
    
@d555 - not in this instance. The original poster had used a sequence of definitions, Xarray[some_Integer] rather than a vector with parts. Both approaches are possible in Mathematica. –  Verbeia Sep 8 '13 at 23:14
    
Ok, sorry. You're right. I get used to use [[..]] all the time for all type of task that I will need and structure indexed. –  d555 Sep 8 '13 at 23:28

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[].

CholeskyDecomposition[] example

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.

share|improve this answer
1  
Jagra, I think you should define sample list, matrix1, tensor1 objects so that this code can be evaluated. I know it's just an illustration, but to someone unfamiliar with what these do, it's not really a useful illustration without that. –  Mr.Wizard Jul 5 '12 at 14:01
    
Noted. Will do. –  Jagra Jul 5 '12 at 14:29
1  
Jagra, I'm sorry that I am criticizing your answer a second time: it seems to me that the second half of your answer is not really appropriate here because any language can have library functions and this says little about the semantics of the language itself or good coding practice. Thanks for updating the first half however, as I think that showing that native arithmetic operations are vector, matrix, and tensor aware is very appropriate. –  Mr.Wizard Jul 6 '12 at 6:46

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.