[R] Avoiding two loops

[Petr Savicky]

On Fri, Feb 04, 2011 at 02:03:22PM -0500, sudhir cr wrote:
> Hello,
> 
> I have a R code for doing convolution of two functions:
> 
> convolveSlow <- function(x, y) {
>     nx <- length(x); ny <- length(y)
>     xy <- numeric(nx + ny - 1)
>     for(i in seq(length = nx)) {
>         xi <- x[[i]]
>         for(j in seq(length = ny)) {
>             ij <- i+j-1
>             xy[[ij]] <- xy[[ij]] + xi * y[[j]]
>         }
>     }
>     xy
> }
> 
> How do I reduce the 2 loops so that I can run the  code faster?

Hello:

Convolution of two vectors may be computed also using matrix reshaping
without a loop. For example 

  convolution <- function(x, y)
  {
      # more efficient if length(x) >= length(y)
      m <- length(x)
      n <- length(y)
      zero <- matrix(0, nrow=n, ncol=n)
      a <- rbind(x %o% y, zero)
      k <- m + n - 1
      b <- matrix(c(a)[1:(n*k)], nrow=k, ncol=n)
      rowSums(b)
  }

Testing this on computing the product of the polynomials (1+t)^4 (1+t)^3
yields

  x <- choose(4, 0:4)
  y <- choose(3, 0:3)
  convolution(x, y)

  [1]  1  7 21 35 35 21  7  1

which is the same as choose(7, 0:7).

See also ?convolve.

Hope this helps.

Petr Savicky.