Quickest way to sum over multiple powers in R-CodePudding

Given two vectors x,y of length d, what is the quickest way to compute

x[1] * y[1]   ...   x[d] * y[d]   x[1]^2 * y[1]^2   ...   x[d]^M * y[d]^M

Both d and M are large and I have a bunch of vectors for which this has to be done

CodePudding user response：

If x and y do not have any components equal to 1 then using the formula for the sum of terms of a geometric series we have the following which reduces the sum from 2M terms to 4 terms. Let z be x*y. Then

z   ...   z^M
= (z^(M 1) - x)/(z - 1)

For example,

M <- 3
x <- y <- 2:5

z <- x * y    
sum( (z^(M 1) - z) / (z - 1) )
## [1] 21546 


# check    
sum(z   z^2   z^3)
## [1] 21546

If x and y can have components equal to 1 then use this:

M <- 3
x <- y <- 1:4

z <- x * y
sum(ifelse(z == 1, M, (z^(M 1) - z)/(z - 1)))
## [1] 5274

sum((z^(M 1) - z   M * (z == 1)) / (z - 1   (z == 1)))
## [1] 5274

# check
sum(x * y   x^2 * y^2   x^3 * y^3)
## [1] 5274

Benchmark

set.seed(123)

x <- runif(10000)
y <- runif(10000)
M <- 100

microbenchmark::microbenchmark(
  A = sum(sapply(seq(M), function(m) x^m * y^m)),
  B = {z <- x * y; sum((z^(M 1)-z M*(z==1))/(z-1 (z == 1))) }
)

Unit: milliseconds
 expr      min       lq       mean    median       uq      max neval cld
    A 209.0281 213.3180 238.843501 221.94970 238.1094 886.1139   100   b
    B   1.3416   1.3703   1.442492   1.39435   1.4391   2.2754   100  a

Update

The original answer used where it should have been *. Have corrected. Note that it now runs even faster.

CodePudding user response：

Speed is unlikely to be a problem here, since you are carrying out simple arithmetic on two vectors. A concise way to write the code would be:

sum(sapply(seq(M), function(m) x^m * y^m))

In terms of speed, this will handle the sum of x and y vectors raised from the first to the hundredth power, even if x and y have a of length 10,000, in about 150 milliseconds:

x <- runif(10000)
y <- runif(10000)
M <- 100

f <- function() {
  sum(sapply(seq(M), function(m) x^m * y^m))
}

microbenchmark::microbenchmark(f())
#> Unit: milliseconds
#>  expr      min       lq     mean   median       uq      max neval
#>   f() 142.1891 145.1826 161.8941 151.1814 163.7553 347.8592   100

^{Created on 2022-04-16 by the reprex package (v2.0.1)}