Skip to content

Inaccuracy in ?convolve

5 messages · Hervé Pagès, Martin Maechler

Original
Hervé Pagès
# 
Hi,

Man page for 'convolve' says:

    conj: logical; if 'TRUE', take the complex _conjugate_ before
          back-transforming (default, and used for usual convolution).

The complex conjugate of 'x', of 'y', of both?

In fact it seems that it takes the complex conjugate of 'y' only which
is OK but might be worth mentioning because (1) conj=TRUE is the default
and (2) with this default then convolve(x,y) is not the same as convolve(y,x).
Note that 'convolve' is commutative with conj=FALSE and would be too
if conj=TRUE was taking the complex conjugate of x _and_ y...

Cheers,
H.
Hervé Pagès
# 
Hi again,

There are many problems with current 'convolve' function.
The author of the man page seems to be aware that 'convolve' does _not_ the
usual thing:

  Note that the usual definition of convolution of two sequences 'x'
  and 'y' is given by 'convolve(x, rev(y), type = "o")'.

and indeed, it does not:

  > x <- 1:3
  > y <- 3:1
  > convolve(x, y, type="o")
  [1]  1  4 10 12  9

The "usual" convolution would rather give:

  > convolve(x, rev(y), type="o")
  [1]  3  8 14  8  3

Also the "usual" convolution is commutative:

  > convolve(y, rev(x), type="o")
  [1]  3  8 14  8  3

but convolve is not:

  > convolve(y, x, type="o")
  [1]  9 12 10  4  1

Of course I could write the following wrapper:

  usual_convolve <- function(x, y, ...) convolve(x, rev(y))

to work around those issues but 'convolve' has other problems:

  (1) The input sequences shouldn't need to have the same length when
      type = "circular" (the shortest can be right-padded with 0s up
      to the length of the longest).
  (2) If the input sequences are both integer vectors, then the result
      should be an integer vector too.
  (3) The "filter" feature seems to be broken (it's not even clear
      what it should do or why we need it though):
        > x <- 1:9
        > y <- 1
        > convolve(x, y, type="f")
        Error in convolve(x, y, type = "f") : subscript out of bounds
        > convolve(y, x, type="f")
        numeric(0)
  (4) If you look at the source code, you'll see that 'x' is first left-padded
      with '0's. The "usual" convolution doesn't do that: it always padd
      sequences on the _same_ side (generally on the right).
  (5) It's not clear why we need a 'conj' arg. All what it does is
      take the conjugate of fft(y) before it does the product with fft(x).
      But this has the "non-usual" effect of reverting the expected result:
        > round(convolve(as.integer(c(0,0,0,1)), 1:7, type="o"))
        [1] 0 0 0 7 6 5 4 3 2 1

Here below is my version of 'convolve' just in case. It does the "usual"
convolution plus:
  - no need to have 'x' and 'y' of the same length when 'type' is "circular",
  - when 'x' and 'y' are integer vectors, the output is an integer vector,
  - no more 'conj' arg (not needed, only leads to confusion),
  - when type is "filter", the output sequence is the same as with
    type="open" but is truncated to the length of 'x' (the original signal)
    It can be seen has the result of 'x' filtered by 'y' (the filter).

convolve2 <- function(x, y, type = c("circular", "open", "filter"))
{
    type <- match.arg(type)
    nx <- length(x)
    ny <- length(y)
    if (type == "circular")
        nz <- max(nx, ny)
    else
        nz <- nx + ny - 1
    if (nz > nx)
        x[(nx+1):nz] <- as.integer(0)
    if (nz > ny)
        y[(ny+1):nz] <- as.integer(0)
    fx <- fft(x)
    fy <- fft(y)
    fz <- fx * fy
    z <- fft(fz, inverse=TRUE) / nz
    if (is.numeric(x) && is.numeric(y))
        z <- Re(z)
    if (is.integer(x) && is.integer(y))
        z <- as.integer(round(z))
    if (type == "filter")
        z[1:nx]
    else
        z
}

Cheers,
H.
Hervé Pagès
# 
Last but not least: convolve2 can be made 100 times or 1000 times faster
than convolve by choosing a power of 2 for the length of the fft-buffer
(a length of 2^n is the best case for the fft, the worst case being when
the length is a prime number):
user  system elapsed
 76.428   0.016  76.445
user  system elapsed
  0.164   0.012   0.179

Here is the modified 'convolve2':

convolve2 <- function(x, y, type = c("circular", "open", "filter"))
{
    type <- match.arg(type)
    nx <- length(x)
    ny <- length(y)
    if (type == "circular") {
        nz <- max(nx, ny)
    } else {
        nz0 <- nx + ny - 1
        nz <- 2^ceiling(log2(nz0))
    }
    if (nz > nx)
        x[(nx+1):nz] <- as.integer(0)
    if (nz > ny)
        y[(ny+1):nz] <- as.integer(0)
    fz <- fft(x) * fft(y)
    z <- fft(fz, inverse=TRUE) / nz
    if (type == "open") {
        z <- z[1:nz0]
    } else {
        if (type == "filter")
            z <- z[1:nx]
    }
    if (is.numeric(x) && is.numeric(y))
        z <- Re(z)
    if (is.integer(x) && is.integer(y))
        z <- as.integer(round(z))
    z
}

In fact, it should try to be smarter than that and not use the fft at all
when one of the 2 input sequences is very short (less than 3 or 4) or
e.g. when one is 10000 times shorter than the other one.

Cheers,
H.
Herve Pages wrote:
Martin Maechler
# 
Thank you, Herve,
Herve> Last but not least: convolve2 can be made 100 times or 1000 times faster
    Herve> than convolve by choosing a power of 2 for the length of the fft-buffer
    Herve> (a length of 2^n is the best case for the fft, the worst case being when
    Herve> the length is a prime number):

    >> x <- 1:100003
    >> y <- 1:1
    >> system.time(cc <- convolve(x, y, type="o")) # uses buffer length of 100003
    Herve> user  system elapsed
    Herve> 76.428   0.016  76.445
    >> system.time(cc <- convolve2(x, y, type="o")) # uses buffer length of 2^17
    Herve> user  system elapsed
    Herve> 0.164   0.012   0.179

The typical approach here and definitely the idea of the
original author of convolve() - would be to use  nextn()  here
instead of  "next_power_of_2()".

convolve() is one of the very old R functions stemming from
Auckland already seen in the very first R tarball that's still
available: Dated from June 20, 1995, with most file dates from
Jun 16/17, i.e. really of even older date,  the file src/rrc/fft
(no 'library', noR extension yet) contains definitions for  'fft', 'nextn'
and 'convolve' where the latter was (note the ":=" for what now
would be assignment to the base package)

convolve := function (a, b, conj=F)
{
	na <- length(a)
	nb <- length(b)
	n <- max(na, nb)
	if (nb < n)
		b <- c(b, rep(0, n - nb))
	else if (na < n)
		a <- c(b, rep(0, n - na))
	da <- fft(a)
	db <- fft(b)
	if(conj) {
		a <- da$x * db$x + da$y * db$y
		b <- - da$x * db$y + db$x * da$y
	}
	else {
		a <- da$x * db$x - da$y * db$y
		b <- da$x * db$y + db$x * da$y
	}
	fft(a, b, inv=T)$x
}

and just for historical fun here's the help file in that
cute olde help file format:

TITLE(convolve @ Fast Convolution)
USAGE(
convolve(a, b, conj=false)
)
ARGUMENTS(
ARG(a,b @ the sequences to be convolved.)
ARG(conj @ logical, if true the transform of LANG(b) is conjugated
before back-transformation.)
)
DESCRIPTION(
LANG(convolve) uses the Fast Fourier Transform to compute the
convolution of the sequences given as its arguments.
PARA
Complex conjugation is useful when computing
autocovariances and autocorrelations by fast convolution.
)
REFERENCES(
Brillinger, D. R. (1981).
ITALIC(Time Series: Data Analysis and Theory), Second Edition.
San Francisco: Holden-Day.
)
SEEALSO(
LANG(LINK(fft)), LANG(LINK(nextn)).
)

Later I had added bits to the docu, convolve got the 'type'
argument, Brian also fixed code and amplified the description and provided
the alternative filter()  which had hencefor been recommended
instead of convolve for most applications.

For back-compatibility (and reverence to the very first R code
writers (!?)) we did not change its behavior but rather documented it
more precisely.

I haven't studied the details of all the possible padding
options for a long time, but I remember that 0-padding
(to length 'n' where 'n' is "highly composite") is only
approximately the same as a non-padded version -- since the
Fourier frequencies  2*pi*j/n depend on n.
As you notice, padding is often very recommendable but it's
strictly giving the solution to a different problem.
For that reason, ?convolve has mentioned  nextn() for 12 years
now, but not "urged" the padding

If we would change convolve() to behave more like your proposal
how many user-defined functions and scripts will give wrong
answers if they are not amended?  
Probably only a very small fraction of all existing code (since
convolve() is not in wide use), but that may still be 100's of
cases. So that would need a the new version with a new name
(such as "convolve2" or maybe slightly nicer "convolve.".

Is it worth introducing that and to start deprecating convolve() ?
IIRC, the last time we considered this (several years ago), we
had concluded "no", but maybe it's worth to get rid of this
infelicity rather than to document it in ever more details.

Martin



    Herve> Here is the modified 'convolve2':

    Herve> convolve2 <- function(x, y, type = c("circular", "open", "filter"))
    Herve> {
    Herve> type <- match.arg(type)
    Herve> nx <- length(x)
    Herve> ny <- length(y)
    Herve> if (type == "circular") {
    Herve> nz <- max(nx, ny)
    Herve> } else {
    Herve> nz0 <- nx + ny - 1
    Herve> nz <- 2^ceiling(log2(nz0))
    Herve> }
    Herve> if (nz > nx)
    Herve> x[(nx+1):nz] <- as.integer(0)
    Herve> if (nz > ny)
    Herve> y[(ny+1):nz] <- as.integer(0)
    Herve> fz <- fft(x) * fft(y)
    Herve> z <- fft(fz, inverse=TRUE) / nz
    Herve> if (type == "open") {
    Herve> z <- z[1:nz0]
    Herve> } else {
    Herve> if (type == "filter")
    Herve> z <- z[1:nx]
    Herve> }
    Herve> if (is.numeric(x) && is.numeric(y))
    Herve> z <- Re(z)
    Herve> if (is.integer(x) && is.integer(y))
    Herve> z <- as.integer(round(z))
    Herve> z
    Herve> }

    Herve> In fact, it should try to be smarter than that and not use the fft at all
    Herve> when one of the 2 input sequences is very short (less than 3 or 4) or
    Herve> e.g. when one is 10000 times shorter than the other one.

    Herve> Cheers,
    Herve> H.
Herve> Herve Pages wrote:
>> Hi again,
    >> 
    >> There are many problems with current 'convolve' function.
    >> The author of the man page seems to be aware that 'convolve' does _not_ the
    >> usual thing:
    >> 
    >> Note that the usual definition of convolution of two sequences 'x'
    >> and 'y' is given by 'convolve(x, rev(y), type = "o")'.
    >> 
    >> and indeed, it does not:
    >> 
    >> > x <- 1:3
    >> > y <- 3:1
    >> > convolve(x, y, type="o")
    >> [1]  1  4 10 12  9
    >> 
    >> The "usual" convolution would rather give:
    >> 
    >> > convolve(x, rev(y), type="o")
    >> [1]  3  8 14  8  3
    >> 
    >> Also the "usual" convolution is commutative:
    >> 
    >> > convolve(y, rev(x), type="o")
    >> [1]  3  8 14  8  3
    >> 
    >> but convolve is not:
    >> 
    >> > convolve(y, x, type="o")
    >> [1]  9 12 10  4  1
    >> 
    >> Of course I could write the following wrapper:
    >> 
    >> usual_convolve <- function(x, y, ...) convolve(x, rev(y))
    >> 
    >> to work around those issues but 'convolve' has other problems:
    >> 
    >> (1) The input sequences shouldn't need to have the same length when
    >> type = "circular" (the shortest can be right-padded with 0s up
    >> to the length of the longest).
    >> (2) If the input sequences are both integer vectors, then the result
    >> should be an integer vector too.
    >> (3) The "filter" feature seems to be broken (it's not even clear
    >> what it should do or why we need it though):
    >> > x <- 1:9
    >> > y <- 1
    >> > convolve(x, y, type="f")
    >> Error in convolve(x, y, type = "f") : subscript out of bounds
    >> > convolve(y, x, type="f")
    >> numeric(0)
    >> (4) If you look at the source code, you'll see that 'x' is first left-padded
    >> with '0's. The "usual" convolution doesn't do that: it always padd
    >> sequences on the _same_ side (generally on the right).
    >> (5) It's not clear why we need a 'conj' arg. All what it does is
    >> take the conjugate of fft(y) before it does the product with fft(x).
    >> But this has the "non-usual" effect of reverting the expected result:
    >> > round(convolve(as.integer(c(0,0,0,1)), 1:7, type="o"))
    >> [1] 0 0 0 7 6 5 4 3 2 1
    >> 
    >> Here below is my version of 'convolve' just in case. It does the "usual"
    >> convolution plus:
    >> - no need to have 'x' and 'y' of the same length when 'type' is "circular",
    >> - when 'x' and 'y' are integer vectors, the output is an integer vector,
    >> - no more 'conj' arg (not needed, only leads to confusion),
    >> - when type is "filter", the output sequence is the same as with
    >> type="open" but is truncated to the length of 'x' (the original signal)
    >> It can be seen has the result of 'x' filtered by 'y' (the filter).
    >> 
    >> convolve2 <- function(x, y, type = c("circular", "open", "filter"))
    >> {
    >> type <- match.arg(type)
    >> nx <- length(x)
    >> ny <- length(y)
    >> if (type == "circular")
    >> nz <- max(nx, ny)
    >> else
    >> nz <- nx + ny - 1
    >> if (nz > nx)
    >> x[(nx+1):nz] <- as.integer(0)
    >> if (nz > ny)
    >> y[(ny+1):nz] <- as.integer(0)
    >> fx <- fft(x)
    >> fy <- fft(y)
    >> fz <- fx * fy
    >> z <- fft(fz, inverse=TRUE) / nz
    >> if (is.numeric(x) && is.numeric(y))
    >> z <- Re(z)
    >> if (is.integer(x) && is.integer(y))
    >> z <- as.integer(round(z))
    >> if (type == "filter")
    >> z[1:nx]
    >> else
    >> z
    >> }
    >> 
    >> Cheers,
    >> H.
    >> 
    >> ______________________________________________
    >> R-devel at r-project.org mailing list
    >> https://stat.ethz.ch/mailman/listinfo/r-devel
    >> 

    Herve> ______________________________________________
    Herve> R-devel at r-project.org mailing list
    Herve> https://stat.ethz.ch/mailman/listinfo/r-devel
1 day later
Hervé Pagès
# 
Hi Martin,

Thanks for taking the time to read me. Here is a followup,
it's rather long but hopefully not too long (and not too boring) ;-)

Quoting Martin Maechler <maechler at stat.math.ethz.ch>:
faster
fft-buffer
being
when
...
The current implementation of convolve uses an fft-buffer of length
nx + ny - 1 for "open" convolution, not nextn().
The fft-based convolution is by nature "circular". However it can be
used for "open" convolutions: the trick is to padd the input sequences
with enough zeros to avoid the overlap inherent to circular convolution.
Then the ouput sequence doesn't "look" circular anymore.

For example, if the input sequences are
   X =   2   5   8
   Y = 100   1   0
then the "circular" convolution of X and Y is
   Z = 205 802 508
To achieve "open" convolution, it's enough to _right_ padd X and
Y with zeros:
   X =   2   5   8   0
   Y = 100   1   0   0
then the "circular" convolution doesn't look circular anymore (it
looks "open"):
   Z = 200 802 508   5

All you need to do is to padd enough zeros to have length X (and Y)
change the result (except that you get extra zeros in the output
sequence):
  X =   2   5   8   0   0   0
  Y = 100   1   0   0   0   0
  Z = 200 802 508   5   0   0

The current approach of adding the strict minimum number of necessary
zeros could _at first glance_ seem reasonable since it minimizes
the amount of memory used. But the problem with this approach is that the
fft algorithm is not very good in the general case and very very bad when
the length of the sequence is a prime number: in this case it has to
perform n^2 complex multiplications so it can litterally take days or
weeks when n is a big prime number (>= 10 millions).
If you ignore the huge bug that 'a' is replaced by 'b' during the 0-padding
operation, this one was better than the current 'convolve'. It already had
the 'conj' argument and the default for it was FALSE so, by default,
the convolve function was doing convolution, not cross-correlation
(more on this below).
But it was doing "circular" convolution only...
filter() might be OK for users who want filtering but what about people
who want a convolution?

Just FYI, with:
   X =   2   5   8
   Y = 100   1
filter(X, Y) gives:
   Z = 502 805  NA
convolve(X, Y, type="filter") gives:
   Z = 208 805
and convolve(x, y, conj=FALSE, type="filter") gives:
   Z = 200 502
so there are some inconsistencies.

But more important: what if I want the "usual" convolution product?

convolve(X, Y, type="o") gives:
   Z = 2 208 805 500
and convolve(x, y, conj=FALSE, type="o") is closer to what I expect
   Z = 5 200 802 508
but still not there :-/

IMO, having a function that does the correct convolution product is
very usefull e.g. for polynom multiplication (if you see X and Y as
the coefficients of polynoms Px and Py, then Z contains the coefficients
of polynomial product Px * Py), or for fast high precision arithmetic
(_exact_ multiplication of numbers with thousands of decimals).
An fft-based convolve function is _the_ tool for those operations
and of course, it shoul be fast.
Nope, it gives exactly the same result (modulo floating point
rounding errors).
convolve() mentions nextn() but does not use it. The user has
no control over the length of the fft-buffer that is used.
But it shouldn't need to: you can safely assume that the best
choice is almost always to take the next_power_of_2 length.
Maybe there could be a few situations where you don't want to do this.
For example, if length(X) + length(Y) - 1 is 1031, the next power
of 2 is 2048 so you need to (almost) double the size of X and Y.
Maybe you are on a machine with very limited memory and you don't
want to do this. But in this case you'd better not try to use the
fft-based convolution at all because when n is a prime number (1031),
then it's better to compute the convolution by just doing:

   Z[i] = sum(k, X[k]*Y[i-k])

You'll do much less multiplications! and they will be on
doubles (or integers) when X and Y are double (or integer) vectors
(with the fft-base convolution, all multiplications are one complexes).

All this to say that it's only worth to use an fft-based convolution
if the length of the fft buffer is a power of 2. And for open
convolution you can _always_ do this. And the cost of it is nothing
compared to the gain in speed (use twice the memory to multiply the
speed by 100, 10000 or more!)

For example, compare the speed of 'filter' (not fft-based) with the
speed of 'convolve2':

  > x <- sample(1:9, 100000, replace=TRUE)
  > y <- rep(1:1, 10000)
  > system.time(z1 <- filter(x, y))
     user  system elapsed
   31.486   0.460  32.478
  > system.time(z2 <- convolve2(x, y, type="o"))
     user  system elapsed
    0.400   0.032   0.454

Note that z1 starts and ends with 4999 NAs. All the values between
(filtered signal?) are also in z2. To be precise:

  all.equal(z1[5000:95000], z2[10000:100000])

is TRUE.
With current version of convolve, they are getting wrong answers
anyway. Nobody on this list wants this

  http://www.sciencemag.org/cgi/content/summary/314/5807/1856

to happen just because of a buggy convolve function in R right?
Here below is another 'convolve2' function (just in case). I've
reintroduced the 'conj' argument after I did some homework
and understood what it was really doing. _What_ it
does is make 'convolve' return the cross-correlation of
X and Y instead of the correlation.
_How_ it does this is by taking the conjugate of fft(y)
before to multiply it by fft(x). Without going into too many
details, I'll just say that this "trick" doesn't play well
with the type="o" option hence the special treatment of
y when 'type' is not "circular" and 'conj' is TRUE.
Also, I've reset the default to FALSE so by default convolve2
does convolution, not cross-correlation.

As for the "filter" type, it seems that the idea of its original
author was to truncate the result of the convolution but it's
not clear to me where this truncation should occur (to the right?
to the left?). IMO, the convolution product of X by Y _is_ the
result of X filtered by Y, it's just that one terminology comes
from maths and the other from signal processing. And there is nothing
to truncate because if the filter (Y) has a length >= 2, then it's
normal to expect a filtered signal longer than X so I would tend to
say that the "filter" feature is a broken concept.

Best,
H.


convolve2 <- function(x, y, conj=FALSE, type=c("circular", "open", "filter"))
{
    type <- match.arg(type)
    nx <- length(x)
    ny <- length(y)
    if (type == "circular") {
        fft_length <- max(nx, ny)
    } else {
        if (conj) {
            y <- rev(y)
            if (is.complex(y))
                y <- Conj(y)
            conj <- FALSE
        }
        nz <- nx + ny - 1
        fft_length <- 2^ceiling(log2(nz))
    }
    if (fft_length > nx)
        x[(nx+1):fft_length] <- as.integer(0)
    if (fft_length > ny)
        y[(ny+1):fft_length] <- as.integer(0)
    fy <- fft(y)
    if (conj)
        fy <- Conj(fy)
    z <- fft(fft(x) * fy, inverse=TRUE) / length(x)
    if (type == "open") {
        z <- z[1:nz]
    } else {
        if (type == "filter")
            z <- z[1:nx]
    }
    if (is.numeric(x) && is.numeric(y))
        z <- Re(z)
    if (is.integer(x) && is.integer(y))
        z <- as.integer(round(z))
    z
}