I am new to Mathematica and I am trying to measure its performance on a fundamentally sequential procedure.

  1. It involves a function STEP that operates on a List and returns an updated List.

  2. Therefore I use a NestList on this operation. The STEP operation involves a sequential calculation per each item in the List. There is no way to avoid this sequential operation without changing the problem.

  3. Each subsequent operation of STEP needs to use the latest List. In fact, it involves a matrix multiplication involving the latest List. I used a table and individual assignments to make this work.

I see no way of rewriting this code, at least easily, but I also see that a naive MATLAB implementation with two FOR loops produces code that is 4-5X faster on my laptop.

Below is a minimal code that shows what I am trying to do, it really slows down for larger NM and stepcount.

Am I missing something obvious?

Edit: The matrix shown below is not identically zero in general, one could think of it is any random matrix whose diagonals are zero. Here, it is a trivial example chosen to show the structure of the code.

NM = 5;
minitial = 2 RandomInteger[{}, NM] - 1.;
Matrix = IdentityMatrix[5] 0;
stepcount = 10^2;
STEP[m_] := 
 Block[{md = m}, 
  Table[md[[i]] = 
    Sign[Tanh[Matrix[[i, All]].md + RandomReal[{-1, 1}]]], {i, NM}]; 
mm = (NestList[STEP[#] &, minitial, stepcount] + 1.)/2 // ArrayPlot

Edit: Just to clarify what I am trying to do, here is how it can be done in MATLAB. Please note that the for loops are unavoidable in this way of thinking.

for ii=1:NT
    for jj=1:NM
            I  = Matrix(jj,:)*m  
            m(jj) = sign (tanh[I]- rand(-1,1))

Maybe there is a more efficient way of doing this in Mathematica than how I implemented it. Hope this clarifies the problem.

  • $\begingroup$ Your Matrix is identically zero, which leads to dramatic simplification if it's true and desired. $\endgroup$
    – Roman
    Commented Jul 7, 2019 at 17:19
  • $\begingroup$ Dear Roman - The Matrix is not identically zero in general, I just added that as an example to show there is matrix multiplication involved. $\endgroup$
    – anon248
    Commented Jul 7, 2019 at 17:33
  • $\begingroup$ I think that memorizing your function STEP[m_]:=STEP[m]= may help in this case since you are calling step with alternating values of +- 1. $\endgroup$
    – cphys
    Commented Jul 8, 2019 at 8:34
  • 1
    $\begingroup$ I have added a standard implementation in MATLAB. Thinking about replacing NestList with Table. $\endgroup$
    – anon248
    Commented Jul 14, 2019 at 18:16
  • 2
    $\begingroup$ I'd say use Compile, since you're dealing with machine float and ints. $\endgroup$
    – Michael E2
    Commented Jul 14, 2019 at 19:15

3 Answers 3


Using Compile is a straightforward way to speed up procedural code based on machine numbers:


NM = 50;
minitial = 2 RandomInteger[{}, NM] - 1.;
Matrix = IdentityMatrix[NM] 0;
stepcount = 10^4;
STEP[m_] := Block[{md = m}, 
   Table[md[[i]] = Sign[Tanh[Matrix[[i, All]].md + RandomReal[{-1, 1}]]],
    {i, NM}]];
mm1 = (NestList[STEP[#] &, minitial, stepcount] + 1.)/2; // AbsoluteTiming
(*  {1.4655, Null}  *)


cf = Compile[{{minitial, _Real, 1}, {Matrix, _Real, 
     2}, {stepcount, _Integer}},
   Block[{md = minitial},
      If[i == 0,  (* i = 0 probably isn't needed *)
       md[[j]],   (* except to conform with NestList *)
       md[[j]] = 
        Sign[Tanh[Matrix[[j, All]].md + RandomReal[{-1, 1}]]]
      {i, 0, stepcount}, {j, Length@minitial}]
    ](*, CompilationTarget -> "C"*)

NM = 50;
minitial = 2 RandomInteger[{}, NM] - 1.;
mm2 = cf[minitial, Matrix, stepcount]; // AbsoluteTiming
(*  {0.162614, Null}  *)

mm1 == mm2
(*  True  *)

Use CompilationTarget -> "C" and it speeds up by another factor of 2.


Here's an improved answer. It takes advantage of the idea of Compilation and your clarification that the vector updates each time a dot product is done. (I should have noticed that). I don't know if it meets your speed requirements, but it might help. I have to use an old form of ReplacePart because the current form won't work in the easy compiler.:(

 nm = 20;
 m = RandomReal[{-1, 1}, {nm, nm}]

 cstep = Compile[{{m, _Real, 2}, {len, _Integer}, {vector, _Real, 
1}, {round, _Integer}},
With[{i = Mod[round, len, 1]}, 
  Sign[Tanh[m[[i]].vector] + RandomReal[{-1, 1}]], i]

 step[m_, len_, {vector_, round_}] := 
 {cstep[m, len, vector, round], round + 1};

 NestList[step[m, nm, #] &, {RandomChoice[{-1, 1}, nm], 1}, 
 steps] // (Map[First] /* (Rescale[#, {-1, 1}, {0, 1}] &) /* ArrayPlot)
  • $\begingroup$ Thank you very much for taking the time -- looks right, I will compare it with the other answer that was given. $\endgroup$
    – anon248
    Commented Jul 15, 2019 at 1:53

Maybe I am missing a subtlety of the problem, but I think there is a relatively simple solution. I can't vouch for speed, but I think the basic issue with your code is that it is needlessly doing element-by-element matrix multiplication when it can all be done at once.

 step = Curry[   
        Function[{matrix, v},   
        Sign[Tanh[matrix.v + RandomReal[{-1, 1}, Length[v]]]]
        {1, 2}];

Then all you need to do is:

 nm = 5;
 m = RandomReal[{-1, 1}, {nm, nm}] (*or whatever you want*);
 minitial = 2 RandomInteger[{}, nm] - 1.;
 NestList[step[m], minitial, 100]//((Rescale[#,{-1,1},{0,1}]&)/*ArrayPlot)

And you get the groovy patterns I assume you are trying to generate.

  • $\begingroup$ It seems the i-th element of md is being changed while the "matrix product" with md is being computed: Table[md[[i]] = Sign[Tanh[Matrix[[i, All]].md.... Isn't that different than the straight Dot[] you're using? Perhaps that's an error on the OP's part, for it seems a strange thing to do. $\endgroup$
    – Michael E2
    Commented Jul 14, 2019 at 17:09
  • 1
    $\begingroup$ @Seth: I haven't completely understood your code but ran it on a single example. Maybe I have not explained the subtlety clearly but my code does not allow the parallel updating of cells. In other words: If we are calculating a quantity Sign[Tanh[matrix.v]+RandomReal[{-1,1}] per each element of "v", this needs to be done sequentially, where every time the matrix calculation is done, the updated list of "v" needs to be used, taking into account the previous calculation. This is needed in what's called "Gibbs Sampling". Your code amounted to updating each "v" in parallel. $\endgroup$
    – anon248
    Commented Jul 14, 2019 at 18:04
  • $\begingroup$ @Michael Perhaps I am carrying MATLAB baggage but Matrix[[i,All]].md just means that one row of the Matrix needs to be used since a single number needs to be calculated to update each m[[i]]. One could calculate the full vector but this would be useless since only ONE element can be updated at a time. This is related to the sequential behavior of the problem. $\endgroup$
    – anon248
    Commented Jul 14, 2019 at 18:10
  • $\begingroup$ @YNSBRYR Just verifying that you didn't want your code to compute Matrix.md, since matrix products are so common. I hadn't heard of Gibbs sampling. $\endgroup$
    – Michael E2
    Commented Jul 14, 2019 at 18:34
  • $\begingroup$ @YNSBRYR -- You are right; I had missed that subtlety of the problem. I think I have another solution though, which I hope to post shortly. $\endgroup$ Commented Jul 15, 2019 at 0:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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