Simple strategy for approx. GP---looking for a reality check.

[Krzysztof Sakrejda]

I recently started using a Stan model that seemed like a good idea but I have not been able to find the right background material to make sure that I'm not off in left field.  I wanted to check if anyone had comments about
this strategy (pasting code didn't work so it's attached and at: https://github.com/sakrejda/gaussianfrankenspline/blob/master/gaussfrankenspline.stan, excuse the name.)

Problem: I want fit a smooth function to some data without needing to worry about the meaning/location of knots, but I have something like 20000 observations so a GP is going to be slow.  Usually the space the GP would need to cover is relatively restricted (e.g.-a year with 365/6 days) so those 20000 observations would be _really_ dense anyway. 

My Stan code says something like:

- Start with the observations y.
- y_i ~ N(mu_i, sigma_obs)
- mu_i is calculated from a spline, with knots p_k
- p_k ~ GP(...)

So roughly speaking this scales as K knots x N obs (as long as K is small enough so the GP part doesn't dominate).  It avoids issues with observations at identical locations, and even my naively coded Stan version works faster than a pure Stan GP for moderate numbers of observations.  As long as there are enough knots it seems like the parameters of the GP are recovered relatively well when I simulate data and fit the true model... and I can still calculate the covariance for any two points.  The questions I have are:

1) Is this sort of model common in some context?  Does it have a name?
2) Any suggestions for making the Stan implementation faster?
3) Is some other approx. GP approach better suited to a Stan implementation currently?

It's a little odd to think of the uncertainty in GP parameters to depend on the number of knots but other than that I haven't found any other issues with this... which makes me uncomfortable.

Thanks for any insight,

Krzysztof

[Andrew Gelman]

Hi, I’m not sure about this one but we’ve been talking about approximate GP’s using inducing points.  I’m cc-ing Louis Mittel who’s working on this.

Louis:  this is another motivation for us to get this project done, then we’d have some code that people can imitate and generalize!

A

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
<gaussfrankenspline.stan>

[Krzysztof Sakrejda]

On Friday, October 31, 2014 10:46:01 PM UTC-4, Andrew Gelman wrote:

Hi, I’m not sure about this one but we’ve been talking about approximate GP’s using inducing points.  I’m cc-ing Louis Mittel who’s working on this.

"inducing points" was the bit of jargon I needed---I think I've found the core lit. now.  Thank you! 
 

Louis:  this is another motivation for us to get this project done, then we’d have some code that people can imitate and generalize!

I'd still love to hear any Stan efficiency comments.  Krzysztof
 

[Bob Carpenter]

Thanks for sharing the model. The code's pretty much optimal
as you wrote it. I also really like the way the names of the functions
and the way they're organized make the model understandable
from the model down.

sq_circ_distance:

rather than making xy a vector, it's a bit faster to use
the two-arg fmin:


if (x > y)
return square(fmin(x - y,
2 * pi() - x + y));
else
return square(fmin(y - x,
2 * pi() - y + x));


The issue is that vectors need to allocate and deallocate memory in
their constructors and destructors and then you have to do the
assignments, then a loop to do the min with more assignments. So doing
something like fmin(xy) with the attendant creation and assignment is
probably 10 to 100 times slower than fmin(x,y). But this is going to
be such a small fraction of your overall time that you probably won't
even notice this difference in the overall program time.

Also square(x) and sqrt(x) are a bit more efficient (and I think
easier to read) than pow(x,2) and pow(x,0.5). The pow() functions
do some conditionals up front and delegate to these if available.
The reason the square() and sqrt() are faster is that they have
precomputed derivatives and don't need to call the expensive pow()
function or store the power. So again, you're unlikely to see any
real speedup.

I like the following even better, because I think it's
clearer what's going on:

real diff;
diff <- fabs(x - y);
return square(fmin(diff, 2 * pi() - diff));

It's maybe a little faster, replacing a binary arithmetic
operation with an fabs() operation, which reduces to either
a negation or a no-op depending on the input sign.

gp_generalized_sq_exp:

I don't know if this is related to the built-in GP density,
which would be more efficient.

model:

you can replace the loop assigning to mu with:

mu <- rep_vector(theta_mu, n_knots);

It's not really going to be any faster, though, just
shorter. But you could then get rid of the local
variable mu altogether and just drop it into multi_normal:

knot_weights ~ multi_normal(rep_vector(theta_mu, n_knots),
theta_Sigma);

Again, you probably won't see any speedup from this, just shorter
code.

If you wrote a function depending on yday_circular_spline that
took a vector day_of_year args, then you could put that loop in
a function and just call it in the final line so it'd look like:

positions ~ normal(yday_circular_spline(day_of_year,
knot_points, knot_weights, knot_scale),
1);

- Bob

[Krzysztof Sakrejda]

On Saturday, November 1, 2014 1:44:09 PM UTC-4, Bob Carpenter wrote:

Thanks for sharing the model.  The code's pretty much optimal
as you wrote it.  I also really like the way the names of the functions
and the way they're organized make the model understandable
from the model down.

You're welcome, thanks for the detailed response.  I don't think there's
really anywhere better to get feedback!  I took your advice and cleaned up
my distance/power calculations.  Being able to write functions is a huge
help for bigger models!

I did run into one issue with overloading
I thought I could just define a new version of yday_circular_spline which
took a vector day_of_year and returned a vector of positions but that wouldn't
link so instead I renamed the non-vectorized version.  Here's the error:

krzysiek@fawkes ~/packages/cmdstan-2.5.0/cmdstan $ make ~/packages/gaussianfrankenspline/gaussfrankenspline

--- Translating Stan model to C++ code ---
bin/stanc /home/krzysiek/packages/gaussianfrankenspline/gaussfrankenspline.stan --o=/home/krzysiek/packages/gaussianfrankenspline/gaussfrankenspline.cpp --no_main
Model name=gaussfrankenspline_model
Input file=/home/krzysiek/packages/gaussianfrankenspline/gaussfrankenspline.stan
Output file=/home/krzysiek/packages/gaussianfrankenspline/gaussfrankenspline.cpp

--- Linking C++ model ---
g++ -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS -I src -I stan/src -isystem stan/lib/eigen_3.2.0 -isystem stan/lib/boost_1.54.0 -Wall -pipe -DEIGEN_NO_DEBUG -Wno-unused-local-typedefs  -lpthread  -O3 -o /home/krzysiek/packages/gaussianfrankenspline/gaussfrankenspline src/cmdstan/main.cpp -include /home/krzysiek/packages/gaussianfrankenspline/gaussfrankenspline.cpp -Lbin -lstan
In file included from <command-line>:0:0:
/home/krzysiek/packages/gaussianfrankenspline/gaussfrankenspline.cpp:457:8: error: redefinition of ‘struct gaussfrankenspline_model_namespace::yday_circular_spline_functor__’
 struct yday_circular_spline_functor__ {
        ^
In file included from <command-line>:0:0:
/home/krzysiek/packages/gaussianfrankenspline/gaussfrankenspline.cpp:422:8: error: previous definition of ‘struct gaussfrankenspline_model_namespace::yday_circular_spline_functor__’
 struct yday_circular_spline_functor__ {
        ^
In file included from stan/src/stan/agrad/rev/operators/operator_unary_plus.hpp:7:0,
                 from stan/src/stan/agrad/rev/operators.hpp:22,
                 from stan/src/stan/agrad/rev.hpp:13,
                 from stan/src/stan/model/model_header.hpp:16,
                 from /home/krzysiek/packages/gaussianfrankenspline/gaussfrankenspline.cpp:3,
                 from <command-line>:0:
stan/src/stan/math/constants.hpp:63:18: warning: ‘stan::math::NEGATIVE_EPSILON’ defined but not used [-Wunused-variable]
     const double NEGATIVE_EPSILON = - std::numeric_limits<double>::epsilon();
                  ^
make/models:9: recipe for target '/home/krzysiek/packages/gaussianfrankenspline/gaussfrankenspline' failed
make: *** [/home/krzysiek/packages/gaussianfrankenspline/gaussfrankenspline] Error 1

Here's the working function block below.  If I remove the _internal suffix it parses but gives the above error during linking.

  real yday_circular_spline_internal(
    real yday,
    vector knot_points, vector knot_weights, real knot_scale
  ) {
    real theta;
    real y;
    theta <- day_to_angle(yday);
    y <- circular_spline(theta, knot_points, knot_weights, knot_scale);
    return y;
  }

  vector yday_circular_spline(
    vector yday,
    vector knot_points, vector knot_weights, real knot_scale
  ) {
    vector[num_elements(yday)] positions;
    for ( i in 1:num_elements(yday)) {
      positions[i] <- yday_circular_spline_internal(
        yday[i], knot_points, knot_weights, knot_scale );
    }
    return(positions);
  }

This was all with Stan 2.5.0 through CmdStan.  Maybe I'm missing something about overloading but I can't tell from the manual.  If it's not clear what the problem is I'll cook up a reproducible example.


 

[Bob Carpenter]

Thanks---that's a known issue

https://github.com/stan-dev/stan/issues/1096

I need to clean up for 2.5.0++. Extending functions to functors allowed
us to plug them into diff eqs, but broke our ability to overload them.

- Bob

[Krzysztof Sakrejda]

On Monday, November 3, 2014 4:12:02 PM UTC-5, Bob Carpenter wrote:

Thanks---that's a known issue

  https://github.com/stan-dev/stan/issues/1096

I need to clean up for 2.5.0++.   Extending functions to functors allowed
us to plug them into diff eqs, but broke our ability to overload them.

- Bob

Ok, thanks for pointing it out.  I just saw the manual and didn't think to look at the issues. Krzysztof

 

[Bob Carpenter]

> On Nov 4, 2014, at 10:07 AM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>
> On Monday, November 3, 2014 4:12:02 PM UTC-5, Bob Carpenter wrote:
> Thanks---that's a known issue
>
> https://github.com/stan-dev/stan/issues/1096
>
> I need to clean up for 2.5.0++. Extending functions to functors allowed
> us to plug them into diff eqs, but broke our ability to overload them.
>
> - Bob
>
> Ok, thanks for pointing it out. I just saw the manual and didn't think to look at the issues. Krzysztof
>

That's our fault --- Stan should do what the manual
says it does! Please feel free to report issues
you see without scouring the issues lists, which are
expansive (because they include refactoring and feature
requests) and hard to search.

- Bob