Vectorization of nested for loop

Question

I have the following function:

g.Bn = function(n) {
  Bn = 0
  for(k in 0:n) {
    res.of.loop = 0
    for(j in 0:k) {
      res.of.loop = res.of.loop + (-1)^j * (j + 1)^n * choose(k, j)
    }
    Bn = res.of.loop * 1/(k+1) + Bn
  }
  return(Bn)
}

Is here a way to vectorize it instead of using for loops?

You could try replacing the inner loop with j<-0:k; res.of.loop<- sum((-1)^j * (j + 1)^n * choose(k, j)) — Dave2e
– Dave2e, Commented Sep 17, 2019 at 20:49
Should modified g.Bn be a vectorized function? Meaning a function that returns a vector of answers when provided with a vector of input? — slava-kohut
– slava-kohut, Commented Sep 17, 2019 at 21:13
@slava-kohut Maybe my question was unclear. I want to use vector operations but the returned value should still be a real number. — JWallmark
– JWallmark, Commented Sep 18, 2019 at 6:25
Oh wait. You could implement in c. That might well lead to substantial speed improvements. — waferthin
– waferthin, Commented Sep 18, 2019 at 7:15

waferthin · Accepted Answer · 2019-09-17 20:59:38Z

4

You could vectorise the inner loop (as per @DaveT), and use sapply:

g.Bn2 = function(n) {
  sum(sapply(0:n, function(k) {
    sum((-1)^(0:k) * (0:k + 1)^n * choose(k, 0:k)) * 1/(k+1)
  }))
}

Or another possibility to vectorise the outer loop:

g.Bn3 = function(n) {
  f <- function(k, n) sum((-1)^(0:k) * (0:k + 1)^n * choose(k, 0:k)) * 1/(k+1)
  sum(Vectorize(f, vectorize.args = "k")(0:n, n))
}

> microbenchmark(g.Bn(100), g.Bn2(100), g.Bn3(100))
       expr      min        lq      mean   median        uq      max neval
  g.Bn(100) 1493.086 1533.9280 1841.3455 1585.354 1675.3575 9023.316   100
 g.Bn2(100)  617.063  650.7850  905.6899  738.230  788.7305 9224.460   100
 g.Bn3(100)  685.094  772.3785 1015.9182  816.945  860.1775 8213.777   100

edited Sep 17, 2019 at 20:59

answered Sep 17, 2019 at 20:52

waferthin

1,6021 gold badge17 silver badges29 bronze badges

Sign up to request clarification or add additional context in comments.

5 Comments

slava-kohut Over a year ago

Tried your code, and g.Bn2 (g.Bn3) produces different results after n passes some threshold. Try g.Bn(15) and g.Bn2(15). Why is that?

waferthin Over a year ago

Honestly, I'm not sure, but I think precision somehow. I can't say which approach is better in this regard.

user20650 Over a year ago

Id guess this may be to do with floating point arithmetic ... try doing it on the log scale?

JWallmark Over a year ago

@waferthin Thanks for replying. I found that just using DaveT's approach for the innerloop and leaving the outer loop as is gives me a better microbenchmark result than any of your solutions. Why do you think that is? Does sapply just iterate in similar fashion to a loop?

g.Bn = function(n) {   Bn = 0   for(k in 0:n) {     j = c(0:k)     res.of.loop = sum((-1)^j * (j + 1)^n * choose(k, j))     Bn = Bn + res.of.loop * 1/(k+1)   }   return(Bn) }

waferthin Over a year ago

Yes the real speed improvement here can be achieved with vectorisation of the inner loop. All other changes are really just different forms of looping (arguably with nicer syntax but no real speed gain). I think the only way to speed up any more would be to think about whether the equations could be simplified in any way.

Yifu Yan · Accepted Answer · 2019-09-17 20:58:24Z

0

You can convert for-loop to map and reduce instead.

In the example below purrr::map iterate all data, and sum reduces a numeric vector into a scaler(a numeric vector with length of one).

g.Bn = function(n) {
    sum(
        purrr::map_dbl(0:n,function(k){
            sum(
                purrr::map_dbl(0:k,function(j){
                    (-1)^j * (j + 1)^n * choose(k, j)
                })
            ) * 1/(k+1)

        })
    )
}

And it seems like all js in inner loop can be replaced with 0:k

g.Bn = function(n) {
    sum(
        purrr::map_dbl(0:n,function(k){
            sum((-1)^(0:k) * (0:k + 1)^n * choose(k, 0:k)) * 1/(k+1)
        })
    )
}

edited Sep 17, 2019 at 20:58

answered Sep 17, 2019 at 20:43

Yifu Yan

6,1561 gold badge21 silver badges37 bronze badges

2 Comments

maydin Over a year ago

Two for replaced by two purrr::. And both are using iteration. What makes it be more powerful? I am not criticizing. I am asking seriously.

Yifu Yan Over a year ago

Using purrr:: saves the trouble of assigning initializer variables at least. And for me, the logic on each level is pretty clear by using nested anonymous function.

Collectives™ on Stack Overflow

Vectorization of nested for loop

2 Answers 2

5 Comments

2 Comments

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

5 Comments

2 Comments

Related