0
$\begingroup$

I would like to perform an iteration on a list, component-wise. The iteration is time consuming and it would be hence beneficial to be able to stop it once an element of the list has reached a threshold after which the computation is unnecessary. To give it some meaning, think for example of a simulation of a many realisations of a Markov chain. Once the absorbing state is reached, it is unnecessary to continue the computstion for thst specific realisation.

A minimum example would be as follows. Having an input list, one adds a list of random integers. The expected result is a list containing the number of iterations for each element to reach a threshold. For example, starting with a state where all the elements equal $-100$ I count the number of iterations for each element to reach $0$.

n = 3;
state = ConstantArray[-100, {n}];
arrivalTime = ConstantArray[0, {n}];
min = 0;
timer = 0;
While[min == 0,  state = state + RandomInteger[{1, 5}, n];
 arrivalTime = 
  With[{\[Lambda] = UnitStep[-state], \[Mu] = Unitize[arrivalTime]}, 
   Subtract[
      1, \[Lambda]] (\[Mu] arrivalTime + 
       Subtract[1, \[Mu]] timer) + \[Lambda] arrivalTime];
 min = Min[arrivalTime];
 ++timer;
 ]

The line

 With[{\[Lambda] = UnitStep[-state], \[Mu] = Unitize[arrivalTime]}, 
   Subtract[
      1, \[Lambda]] (\[Mu] arrivalTime + 
       Subtract[1, \[Mu]] timer) + \[Lambda] arrivalTime];

is a construct I learned on this site, courtesy of HenrikSchumacher. It implements the logic, IF the current state for the element is >=0 and the arrival time is still $0$, as initialised, then write the number of iterations.

This works, but it has the drawback of the update, in this case to compute and add the random integer, to be still performed. It would be great if this could be avoided for elements that have already reached the threshold, in cases where the update is time expensive. One could write for example a for loop over the elements, place an IF clause and proceed with the update only if required, but from what I understand this is not optimal.

What is the best way to achieve this purpose in Mathematica? Thanks a lot

EDIT

Following from the Daniel Huber's answer, I attach a timing comparison. Using 1000 long lists, first the approach in the answer, then what is in the main post.

AbsoluteTiming[inp = Table[-1000, {1000}];
 thresh = 0;
 (sum = #; count = 0;
    While[sum < thresh, count++; sum += RandomInteger[{0, 2}]];
    count) & /@ inp]
(*1.85, {results} *)
AbsoluteTiming[n = 1000;
 state = ConstantArray[-1000, {n}];
 arrivalTime = ConstantArray[0, {n}];
 min = 0;
 timer = 0;
 While[min == 0,  state = state + RandomInteger[{0, 2}, n];
  arrivalTime = 
   With[{\[Lambda] = UnitStep[-state], \[Mu] = Unitize[arrivalTime]}, 
    Subtract[
       1, \[Lambda]] (\[Mu] arrivalTime + 
        Subtract[1, \[Mu]] timer) + \[Lambda] arrivalTime];
  min = Min[arrivalTime];
  ++timer;
  ]
 ]
(* 0.07 , *)

So the approach I am using seems to be much faster, even if the update is performed on every element regardless of having reached the threshold. I am wondering how fast it could be made, if the iterations were stopped for elements that have already reached the threshold. Daniel Huber's code does that, but seems much slower overall, for reasons that I would like to fully understand, thanks a lot.

$\endgroup$
4
  • $\begingroup$ One thing to keep in mind that at least to a certain degree, the speedup of the UnitStep-version of If comes from the fact that it is branchless - this necessarily means that both branches need to be computed every time. (Which is of course part of the tradeoff - you have to check whether the slowdown due to the extraneous computations or the slowdown due to the branching is worse in each case) $\endgroup$
    – Lukas Lang
    Commented Feb 7, 2022 at 22:18
  • $\begingroup$ @Lukas Lang, so there is no way to get the best of both worlds? I need to specifically evaluate if the branching cost is worth the lower number of computations, no other way? $\endgroup$
    – Smerdjakov
    Commented Feb 7, 2022 at 23:21
  • $\begingroup$ @Smerdjakov if you already have an iterative version, use Compile and you get the speedup for free $\endgroup$
    – b3m2a1
    Commented Feb 8, 2022 at 0:29
  • 1
    $\begingroup$ @Smerdjakov in principle yes. However, the good news is that in Mathematica, there are usually several performance wins to be had for way less effort. This includes clever use of auto-vectorizing functions (such as your UnitStep version) and Complie (as suggested by b3m2a1) $\endgroup$
    – Lukas Lang
    Commented Feb 8, 2022 at 7:44

1 Answer 1

1
$\begingroup$

I hope I understand correctly what you want: Given and input vector, for every element count how many time an random integer is added until a given threshold is reached.

This can be implemented by: The input:

inp = Table[-100, {10}];

Then we define a function, that adds random integers to the input until the threshold is reached and returns a count of the additions. Then we apply this function to every element of the input vector:

thresh = 0;
(sum = #; count = 0; 
   While[sum < thresh, count++; sum += RandomInteger[{0, 9}]]; 
   count) & /@ inp
(* {22, 25, 26, 23, 26, 25, 18, 23, 23, 19} *)
$\endgroup$
4
  • $\begingroup$ thanks lot for this. I thought about this approach but it seems it is rather slow, I added a timing comparison in the post, would really like to understand how it works. Does a function like yours get computed on each element in turn? $\endgroup$
    – Smerdjakov
    Commented Feb 7, 2022 at 22:07
  • $\begingroup$ Map (for short: /@) on a vector applies the function to each element. $\endgroup$ Commented Feb 8, 2022 at 9:01
  • $\begingroup$ thanks, I understand what map does in terms of output. I was just wondering if, "under the bonnet", it would use vectorised operations, all in one sweep. Thanks $\endgroup$
    – Smerdjakov
    Commented Feb 8, 2022 at 13:16
  • $\begingroup$ As long as you are not using Parallel... it is always done sequentially. However, with Map, interpretation is only done once. $\endgroup$ Commented Feb 8, 2022 at 14:10

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.