------------------------------
Message: 8
Date: Sun, 4 Feb 2007 17:51:22 -0800
From: hpages at fhcrc.org
Subject: Re: [Rd] convolve: request for "usual" behaviour + some
improvements + some fixes
To: Martin Maechler <maechler at stat.math.ethz.ch>
Cc: r-devel at r-project.org
Message-ID: <1170640282.45c68d9aa3a30 at webmail.fhcrc.org>
Content-Type: text/plain; charset=ISO-8859-1
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>:
"Herve" == Herve Pages <hpages at fhcrc.org>
on Fri, 02 Feb 2007 21:30:04 -0800 writes:
Herve> Last but not least: convolve2 can be made 100 times or 1000 times
Herve> than convolve by choosing a power of 2 for the length of the
Herve> (a length of 2^n is the best case for the fft, the worst case
Herve> the length is a prime number):
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()".
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)
= nx + ny - 1. But you can add more zeros if you want: this doesn't
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).
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)).
)
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...
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.
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.
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.
Nope, it gives exactly the same result (modulo floating point
rounding errors).
For that reason, ?convolve has mentioned nextn() for 12 years
now, but not "urged" the padding
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.
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.
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
}