# Speeding up an intrinsically sequential operation

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.

ClearAll["Global*"]
SeedRandom;
NM = 5;
minitial = 2 RandomInteger[{}, NM] - 1.;
Matrix = IdentityMatrix 0;
stepcount = 10^2;
STEP[m_] :=
Block[{md = m},
Table[md[[i]] =
Sign[Tanh[Matrix[[i, All]].md + RandomReal[{-1, 1}]]], {i, NM}];
md]
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))
end
mm(:,ii)=m;
end


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

• Your Matrix is identically zero, which leads to dramatic simplification if it's true and desired. – Roman Jul 7 '19 at 17:19
• Dear Roman - The Matrix is not identically zero in general, I just added that as an example to show there is matrix multiplication involved. – anon248 Jul 7 '19 at 17:33
• 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. – cphys Jul 8 '19 at 8:34
• I have added a standard implementation in MATLAB. Thinking about replacing NestList with Table. – anon248 Jul 14 '19 at 18:16
• I'd say use Compile, since you're dealing with machine float and ints. – Michael E2 Jul 14 '19 at 19:15

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

OP's:

SeedRandom;
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}  *)


Compiled:

cf = Compile[{{minitial, _Real, 1}, {Matrix, _Real,
2}, {stepcount, _Integer}},
Block[{md = minitial},
Rescale@Table[
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"*)
];

SeedRandom;
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;
steps=1000;
m = RandomReal[{-1, 1}, {nm, nm}]

cstep = Compile[{{m, _Real, 2}, {len, _Integer}, {vector, _Real,
1}, {round, _Integer}},
With[{i = Mod[round, len, 1]},
ReplacePart[vector,
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)

• Thank you very much for taking the time -- looks right, I will compare it with the other answer that was given. – anon248 Jul 15 '19 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.

• 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. – Michael E2 Jul 14 '19 at 17:09
• @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. – anon248 Jul 14 '19 at 18:04
• @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. – anon248 Jul 14 '19 at 18:10
• @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. – Michael E2 Jul 14 '19 at 18:34
• @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. – Seth Chandler Jul 15 '19 at 0:01