ESS calcs?

[Bob Carpenter]

Is anyone working on

a) speeding up ESS calcs and
b) making sure they work for single chains?

Last I heard, that was hugely slowing down
RStan. I'm also not sure if the basic I/O for
the chains themselves is super slow, in which case
that will need to be dealt with.

So before RStan is usable, we either have to

a) fix the C++ or
b) rewrite Andrew and Kenny's ESS and maybe
all the transpositions and index-munging in R.

We saw from profiling the bottleneck is in
the index munging in I/O and in the autocorrelation
functions when inputs are the "wrong" size. For the
sizing, Matt thought we could pad, but might have
to adjust the outputs.

I'm just writing to make sure this doesn't fall
between the cracks. I can look into it when I get
back if no one else wants to.

- Bob

[Jiqiang Guo]

I did not work on this. But I did try to understand a bit of the
business behind
the calculation of auto-correlation using FFT.

My thinking is one options is that we use FFTW in RStan if possible,
for which I do not think
we have a license issue. One reason is that part of the
mcmc::chains's is in RStan; I am
not using mcmc::chains directly any more.

If we use kissFFT and pad more zeros to make the size be good for it,
we need to figure out
how many zeroes we need to add. The transforming of output should be easy.

I did not understand the index munging in I/O mentioned below.

--
Jiqiang

[Bob Carpenter]

On 7/30/12 2:03 PM, Jiqiang Guo wrote:
> I did not work on this. But I did try to understand a bit of the
> business behind
> the calculation of auto-correlation using FFT.
>
> My thinking is one options is that we use FFTW in RStan if possible,

We can just use R at that point.

> for which I do not think
> we have a license issue.

Correct.

> One reason is that part of the
> mcmc::chains's is in RStan; I am
> not using mcmc::chains directly any more.

Are you extending mcmc::chains or did you just
rewrite it? We should talk about the design, because
the whole point of writing mcmc::chains was so that
you could use it in RStan.

> If we use kissFFT and pad more zeros to make the size be good for it,
> we need to figure out
> how many zeroes we need to add.

Powers of 2 work well, so that's no problem. The issue
is whether that changes the output in some predictable
way we can rescale automatically.

An alternative is to do the direct calculations, which
can be expensive if the autocorrelation time is long.

> The transforming of output should be easy.
>
> I did not understand the index munging in I/O mentioned below.

The problem is that C++ is row-major (first-index major for
larger numbers of dimensions) and R is column-major (last-index
major), so you have to write logic to do this transformation.
Right now, that's taking about half the overall time of the
ESS calcs.

- Bob

[Marcus Brubaker]

You can pad out the signal, do the FFT, multiply with the conjugate signal, invert the FFT and then truncate back down to the original signal length.  Since the length of the signal should be much longer than the correlation length of the process, you'll never end up looking at correlations that might be influenced by the zero padding.

I've attached a bit of matlab code to do this which was sent to me by Radford years back. Given the advances with FFTW over the years, it's probably not necessary with Matlab anymore, but it shows how to do it, including getting the right length.

Cheers,
Marcus

[Jiqiang Guo]

On Mon, Jul 30, 2012 at 3:27 PM, Bob Carpenter <ca...@alias-i.com> wrote:

> Are you extending mcmc::chains or did you just
> rewrite it? We should talk about the design, because
> the whole point of writing mcmc::chains was so that
> you could use it in RStan.

I guess we have talked about this before. Essentially I adapted the
split Rhat and effective samples size calculation from mcmc::chains to
two functions implemented in C++, which are in binary form in RStan.
As we talked before, I am trying to minimize the usage of Rcpp's
functionality. So I do not like that a chunk of data sitting in a C++
class, on which I do not have direct control. Previously, I have an
mcmc::chains instance, which I can only control through Rcpp's module
mechanism. Using this method, I could not (or it is difficult for me
to) save the samples and reload them such that all the functions for
mcmc::chains are still usable.

Now all the samples are written directly as objects of R, on which I
could call those two functions to calculated Rhats and ESS regardless
of whether it is another session of R or not. I know this is opposite
to the idea of using mcmc::chains. But in my opinion, when using C++
with R, it is better to just have methods implemented in C++. And we
can talk about this; maybe we have better solutions.

--
Jiqiang

[Bob Carpenter]

I knew there was a problem with serializing the C++ data
structures. I never looked into how to do that directly.
I find the R and Rcpp docs pretty much impenetrable.

If you represent the data from the chains in R, we can
just get rid of Rcpp altogether. We don't want to have
to represent the data in both C++ and R data structures.

It'll be more scalable to write the data to a file
and then read it into R rather than reading it into a C++
data structure then copying it into R (you could delete the
C++ data structure after reading into R, but the peak memory
requirement will be the same unless we can read and delete
incrementally).

You could probably use an FFT from R wrapping FFTW, too.

At that point, why not just go with the very simple
wrapper that I wrote? The main motivation for Rcpp was
to keep all the data and calcs on the C++ side for speed
and scalability.

Then we can write the same kind of functionality again
in C++ for command-line use and not be constrained by
interacting with R (like trying to swap row-major vs.
column-major representations and avoiding 64-bit integer
types).

- Bob

[Jiqiang Guo]

I have patched the function of computing the auto-correlations
by padding more zeros before calling kissFFT.

Now it is a lot faster for 52012, which happens to be 13003 * 4 and 13003
is a prime number.

--
Jiqiang

On Mon, Jul 30, 2012 at 1:50 AM, Bob Carpenter <ca...@alias-i.com> wrote:

[Bob Carpenter]

On 8/2/12 11:17 PM, Jiqiang Guo wrote:
> I have patched the function of computing the auto-correlations
> by padding more zeros before calling kissFFT.

Great news. I hope that will be fast enough.

Any idea how much slower it is than FFTW? (Is there
an FFTW wrapper in R?) The reports I've seen say about
a factor of 2 or 3.

> Now it is a lot faster for 52012, which happens to be 13003 * 4 and 13003
> is a prime number.

In trying to understand why Jiqiang thought to report this,
I just took a look at the code in:

stan/lib/eigen_3.1.0/unsupported/Eigen/src/FFT/ei_kissfft_impl.h

and there are several places where the code really wants the
input size to be factorizable into powers of 2, 3 and 5.

I also looked at the current KissFFT code, and it's changed significantly
since Eigen imported a version (in what looks to be 2009).
So, if the latest patch isn't fast enough, we can try running
the latest KissFFT for comparison.

We'd still have to pad, because it looks like it still wants
the size to be factorizable into powers of 2, 3 and 5. (I have
no idea why -- presumably you could keep going, but it'd take
an awful lot of unfolding until we got to powers of 13003!)

- Bob

[Jiqiang Guo]

Please see below for not that important comments. 

Any idea how much slower it is than FFTW?  (Is there
an FFTW wrapper in R?)  The reports I've seen say about
a factor of 2 or 3.

Using the Eigen wrapped kissFFT or FFTW, the same code (but compiled/linked differently) using kissFFT takes about twice as much time as using FFTW, both compiled with -O3.  But even when using Eigen wrapped FFTW, it is also better to pad more zeros to have a good size though not sure what would happen in other forms of using FFTW.  

There is an fftw R package, but I did not try. 

Now it is a lot faster for 52012, which happens to be 13003 * 4 and 13003
is a prime number.

In trying to understand why Jiqiang thought to report this,
I just took a look at the code in:

I just thought it is a bit funny that we have a large prime number here; otherwise, we might not notice this issue. 

No other thoughts.