Re: [stan-dev/math] Matrix exponential (#32)

[Bob Carpenter]

[off issue --- we want these general discussions on stan-dev,
not on issues]

Once you have a branch, you don't need to outline your changes.
We can just go to GitHub and compare to develop to see what
you've changed.

Rather than the kind of quick checks you're doing, you need
to learn how to write Stan unit tests. Once you do that, they're
even easier, because all the library paths are included.

You can't include the model_header.hpp file, as it's not
part of stan::math. Just try to include only the headers
that you need.

You almost certainly want to be more specific with those
matrix types. Eigen itself uses the curiously recursive
template pattern (CRTP), but we don't need to go that far.

There shouldn't be diffs in the make files. You should
pull the most recent version of develop and make sure you're
in synch. You want to be careful not to unintentionally check
in files---it takes some getting used to.

> stan::math::var E11=-1, E22=-17, lp11;
>
> E_diag << E11, 0,
> 0, E22;

You can just use E_diag directly --- you don't need the intermediate
stan::math::var. and once you're in the stan::math namespace (where
you want to define things you're doing), you don't need the stan::math::
qualifier.

Don't commit files that only have comments, like mat.hpp as is.
We don't want comments left around in our code.

Minor point, but please put space around operators. Check out
the Google C++ style guide, which we're trying to follow.

The doc for frexp mentions sine and I have no idea what's up with
that. Please don't just copy things blindly---this is a common
bad practice and one of the reasons I don't trust most doc---too many
cut and paste errors. Also, you want to reference the <cmath> version,
which should be std::frexp, not ::frexp (which uses the <math.h> version
in the top-level namespace).

You want to declare short functions like matrix_exp as inline.

You don't need Scalar(2) in pow for Stan---it's more efficient
to just send a doulbe into stan::math::pow as the first argument
as it computes fewer derivatives.

Now I can see the problem with their class for your application.
They declared everything as private rather than protected, not assuming
anyone would extend the damn thing. One option might be to get
them to change their class going forward. I don't know if it's
still being used.

Otherwise, we don't want to actually edit their class, you want
to fork it and fix it all in stan::math to be the way you want.
Their library lets you do that---you just have to keep their
copyright and mention at the top of the file that you edited it.

Hope that's not too much feedback from the get go!

- Bob

> On Aug 18, 2016, at 4:32 PM, Charles Margossian <notifi...@github.com> wrote:
>
> @syclik: Yes, I have permission to push. GitHub gymnastic is a little new to me, so I had to set up a few things (following the guidelines here: https://github.com/stan-dev/stan/wiki/Dev:-Git-Process). Sorry for the delay is responding!
>
> Let me start by sharing the non-ideal working prototype, where the exp function works for rev AD (this is easily extendable to fwd AD), and where I hacked into Eigen/unsupported.
>
> Created branch: feature/issue-32-matrix-exponential
>
> Changes from dev:
>
> modified: lib/eigen_3.2.8/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
> modified: make/libraries
> modified: makefile
> modified: stan/math/rev/mat.hpp
> new file: stan/math/rev/mat/fun/matrix_exp.hpp
> modified: stan/math/rev/scal.hpp
> new file: stan/math/rev/scal/fun/frexp.hpp
>
> When I tried to run A.exp(), where A is a matrix of var, there were two functions that did not recognize stan variables: frexp and computeUV. These are the two functions I overloaded. Frexp was from cmath, so this was easy to do without changing Eigen/unsupported. ComputeUV was a class function from MatrixExponential, so I changed it in the class. Subsequently, exp() worked on matrices of var.
>
> I'm still not sure what you mean by "act on private members".
> Do you mean private member variables (non-static variables declared
> in the class) or non-static class functions? If it's
> member variables, are they things you pass into the class to
> start with or things that are easy to recreate?
>
> I mean private member variables. They are passed (in an initialization list) into the class, but they are easy to recreate, seeing they correspond to various properties of the input matrix. My plan is to write a function that creates these variables, and then applies various function members of the MatrixExponential class. Which is almost the same as rewriting the class, and should do the job. What I was looking for was a more elegant solution that took advantage of the MatrixExponential class.
>
> Here's the C++ code I used to do some quick checks on the prototype I posted. In this example, the exp function does not have the intended format, namely exp(A), rather A.exp().
>
> #include <unsupported/Eigen/MatrixFunctions>
> #include <stan/math.hpp>
>
> int main() {
>
>
> // Look at simple of a 1 x 1 matrix
> stan::math::var mu = 2.0, lp;
> Eigen::Matrix<stan::math::var, 1, 1> mu_matrix;
> Eigen::Matrix<stan::math::var, 1, 1> mu_exp;
>
> mu_matrix << 2*mu;
> mu_exp = mu_matrix.exp();
>
> lp = mu_exp(0);
> std::cout << "mu = " << mu.val() << "\n";
> std::cout << "exp(2*mu) = " << lp.val() << "\n";
> lp.grad();
> std::cout << "d.f / d.mu = " << mu.adj() << "\n\n";
>
> // Expected result: d.f / d.mu = 2*exp(4) ~ 110
>
> std::cout << "-------------------------------------------- \n\n";
>
>
> // Look at non-diaginal case (example from Moler & Van Loan, 2003)
> Eigen::Matrix<stan::math::var, 2, 2> E, E_diag, E_b, E_b_1;
> stan::math::var E11=-1, E22=-17, lp11;
>
> E_diag << E11, 0,
> 0, E22;
>
> E_b << 1, 3,
> 2, 4;
>
> E_b_1 << -2, 1.5,
> 1, -.5;
>
> E = E_b*E_diag*E_b_1;
>
> std::cout << "The Matrix E is:\n" << E << "\n\n"
> << "The matrix exponential of E is:\n" << E.exp() << "\n\n";
>
> Eigen::Matrix<stan::math::var, 2, 2> expE = E.exp();
> lp11 = expE(1,1);
> lp11.grad();
> std::cout << "d.exp(E)[2,2]/dx = " << E11.adj() << "\n\n";
> // Expected Result (worked out by hand) for d.exp(E)/dx:
> //
> // dE/dx = -.736 .552
> // -1.472 1.104
> //
> // Note: derivative is identical to E.
>
> return 0;
>
> }
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub, or mute the thread.
>

[Ben Goodrich]

On Thursday, August 18, 2016 at 11:44:12 AM UTC-4, Bob Carpenter wrote:

You can't include the model_header.hpp file, as it's not
part of stan::math.  Just try to include only the headers
that you need.

Also, do not modify Stan's local copy of Eigen or Boost unless absolutely necessary because RStan doesn't use them. Just copy Eigen's unsupported matrix exponential file into stan/prim/mat and then tweak it.

Ben

[Charles Margossian]

They declared everything as private rather than protected, not assuming
anyone would extend the damn thing.

yeah, I too find this upsetting. 

Otherwise, we don't want to actually edit their class, you want
to fork it and fix it all in stan::math to be the way you want.
Their library lets you do that---you just have to keep their
copyright and mention at the top of the file that you edited it.

I think this should work. My only concern is that there are some dependencies in supported Eigen. For example the function MatrixBase<Derived>::exp() is forward declared in MatrixBase.h (line 455) -- but I think there is a way around it.

Hope that's not too much feedback from the get go!

Not at all, thank you for taking the time to write it all out! I'll be sure to put it all in practice.

Also, do not modify Stan's local copy of Eigen or Boost unless absolutely necessary because RStan doesn't use them. Just copy Eigen's unsupported matrix exponential file into stan/prim/mat and then tweak it.

That's the plan!  

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

[Bob Carpenter]

Are you on the stan-dev list? If not, you should be. I'm
cc-ing you just in case.

> On Aug 18, 2016, at 9:19 PM, Charles Margossian <charles.ma...@gmail.com> wrote:
>
> They declared everything as private rather than protected, not assuming
> anyone would extend the damn thing.
>
> yeah, I too find this upsetting.

File an issue with Eigen. They've already fixed a bunch of
stuff we've asked about, as has Boost. That's the beauty of
open source. Even better, submit a pull request for them (I have
no idea what version control system they use).

> Otherwise, we don't want to actually edit their class, you want
> to fork it and fix it all in stan::math to be the way you want.
> Their library lets you do that---you just have to keep their
> copyright and mention at the top of the file that you edited it.
>
> I think this should work. My only concern is that there are some dependencies in supported Eigen. For example the function MatrixBase<Derived>::exp() is forward declared in MatrixBase.h (line 455) -- but I think there is a way around it.

We'll still have the Eigen library to include, so that
won't be a problem.

>
> Hope that's not too much feedback from the get go!
>
> Not at all, thank you for taking the time to write it all out! I'll be sure to put it all in practice.

I type quickly, which is a blessing and a curse. I used to be a
secretary and then a keypunch operator back in the late 70s and
early 80s to earn pocket money.

> Also, do not modify Stan's local copy of Eigen or Boost unless absolutely necessary because RStan doesn't use them. Just copy Eigen's unsupported matrix exponential file into stan/prim/mat and then tweak it.
>
> That's the plan!

But as Ben points out, it's OK to fork them into the
stan::math namespace with attribution if we can't get
them to change Boost.

- Bob

[Charles Margossian]

File an issue with Eigen.  They've already fixed a bunch of
stuff we've asked about, as has Boost.  That's the beauty of
open source.  Even better, submit a pull request for them (I have
no idea what version control system they use).

They use GitHub: https://github.com/RLovelett/eigen

Their latest stable release is 3.2.9, and Stan currently uses 3.2.8. The landscape in MatrixExponential.h looks quite different! Most importantly, not everything is private. I'll try to use the new version of code. 

Is it reasonable to expect Stan will use version 3.2.9 of Eigen? 

On Thu, Aug 18, 2016 at 4:39 PM, Bob Carpenter <ca...@alias-i.com> wrote:

Are you on the stan-dev list?  If not, you should be.  I'm
cc-ing you just in case.

[Bob Carpenter]

Stan is unfortunately pegged to RcppEigen from R, so whenever
they update, we'll update. We're ready for 3.3.0, which will
hopefully be released soon.

Looks like they use Mecurial (that version control system lost
the battle---nobody supports it any more), because that repo
says at the top:

Git mirror of upstream Mercurial based Eigen repository.

By "upstream", I think they mean changes to Eigen's Mercurial
repo get mirrored on GitHub, but probably not vice-versa. You
can see there aren't any pull requests or issues here.

- Bob

> On Aug 19, 2016, at 1:40 PM, Charles Margossian <charles.ma...@gmail.com> wrote:
>
> File an issue with Eigen. They've already fixed a bunch of
> stuff we've asked about, as has Boost. That's the beauty of
> open source. Even better, submit a pull request for them (I have
> no idea what version control system they use).
>
> They use GitHub: https://github.com/RLovelett/eigen
>
> Their latest stable release is 3.2.9, and Stan currently uses 3.2.8. The landscape in MatrixExponential.h looks quite different! Most importantly, not everything is private. I'll try to use the new version of code.
>
> Is it reasonable to expect Stan will use version 3.2.9 of Eigen?
>
>
>
> On Thu, Aug 18, 2016 at 4:39 PM, Bob Carpenter <ca...@alias-i.com> wrote:
> Are you on the stan-dev list? If not, you should be. I'm
> cc-ing you just in case.
>

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

> For more options, visit https://groups.google.com/d/optout.
>
>

> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

[Bob Carpenter]

I'm not sure if RcppEigen contains the unsupported modules.
If not, we can just copy them into Stan so they'll work in
all the interfaces.

- Bob

[Sebastian Weber]

If you ask Dirk, the maintainer of RcppEigen, he is usually quite open for updates. Even better you provide a good reason to update.

Sebastian

On Friday, August 19, 2016 at 2:13:12 PM UTC+2, Bob Carpenter wrote:
> I'm not sure if RcppEigen contains the unsupported modules.
> If not, we can just copy them into Stan so they'll work in
> all the interfaces.
>
> - Bob
>
> > On Aug 19, 2016, at 2:08 PM, Bob Carpenter
> >

> > Stan is unfortunately pegged to RcppEigen from R, so whenever
> > they update, we'll update. We're ready for 3.3.0, which will
> > hopefully be released soon.
> >
> > Looks like they use Mecurial (that version control system lost
> > the battle---nobody supports it any more), because that repo
> > says at the top:
> >
> > Git mirror of upstream Mercurial based Eigen repository.
> >
> > By "upstream", I think they mean changes to Eigen's Mercurial
> > repo get mirrored on GitHub, but probably not vice-versa. You
> > can see there aren't any pull requests or issues here.
> >
> > - Bob
> >
> >> On Aug 19, 2016, at 1:40 PM, Charles Margossian
> >>

> >> File an issue with Eigen. They've already fixed a bunch of
> >> stuff we've asked about, as has Boost. That's the beauty of
> >> open source. Even better, submit a pull request for them (I have
> >> no idea what version control system they use).
> >>
> >> They use GitHub: https://github.com/RLovelett/eigen
> >>
> >> Their latest stable release is 3.2.9, and Stan currently uses 3.2.8. The landscape in MatrixExponential.h looks quite different! Most importantly, not everything is private. I'll try to use the new version of code.
> >>
> >> Is it reasonable to expect Stan will use version 3.2.9 of Eigen?
> >>
> >>
> >>
> >> On Thu, Aug 18, 2016 at 4:39 PM, Bob Carpenter

> >> Are you on the stan-dev list? If not, you should be. I'm
> >> cc-ing you just in case.
> >>

> >>> On Aug 18, 2016, at 9:19 PM, Charles Margossian rote:

[Bob Carpenter]

I followed your lead and submitted the request
as an issue:

https://github.com/RcppCore/RcppEigen/issues/36

explaining we needed the updates to matrix exponentials.

Thanks for the tip.

- Bob

[Joshua N Pritikin]

On Fri, Aug 19, 2016 at 02:46:55PM +0200, Bob Carpenter wrote:
> I followed your lead and submitted the request
> as an issue:
>
> https://github.com/RcppCore/RcppEigen/issues/36
>
> explaining we needed the updates to matrix exponentials.

There is no real urgency to get the new version in RcppEigen. You can
copy the source code into stan until Eigen/RcppEigen catch up.

[Bob Carpenter]

That is indeed the backup plan. It's just more work for
us in the math library and the three interfaces we support,
all of which deal with the libraries slightly differently.

- Bob

[Ben Goodrich]

On Friday, August 19, 2016 at 8:13:12 AM UTC-4, Bob Carpenter wrote:

I'm not sure if RcppEigen contains the unsupported modules.

Yes.

[Daniel Lee]

Yes, it is contained?

[Bob Carpenter]

Thanks. I checked before submitting the RcppEigen issue to
request an update to the latest stable version and should've
posted the answer here.

- Bob

[Ben Goodrich]

On Friday, August 19, 2016 at 11:48:20 AM UTC-4, Daniel Lee wrote:

Yes, it is contained?

Yes

[Charles Margossian]

Stan is unfortunately pegged to RcppEigen from R, so whenever
they update, we'll update.  We're ready for 3.3.0, which will
hopefully be released soon.

Version 3.3.0 of what? 

Latest stable release of Eigen is 3.2.9

Latest release of Rcpp Eigen 0.3.2.8.1

 

--

You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

[Ben Goodrich]

On Friday, August 19, 2016 at 1:15:24 PM UTC-4, Charles Margossian wrote:

Stan is unfortunately pegged to RcppEigen from R, so whenever
they update, we'll update.  We're ready for 3.3.0, which will
hopefully be released soon.

Version 3.3.0 of what? 

Eigen

[Daniel Lee]

See here: http://eigen.tuxfamily.org/index.php?title=Main_Page

[Rob Trangucci]

3.3 looks to be a major update to Eigen according to here: http://eigen.tuxfamily.org/index.php?title=3.3

Rob

[Daniel Lee]

It is.

[Daniel Lee]

Bob and Ben have been working through the math library to get it to compile.

[Ben Goodrich]

On Friday, August 19, 2016 at 4:42:55 PM UTC-4, Daniel Lee wrote:

Bob and Ben have been working through the math library to get it to compile.

Bob has been working. I have been harping.

[Bob Carpenter]

It's all done, right? And isn't it Daniel that figured
out all the traits to make it work? Not that we're
keeping score.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

[Charles Margossian]

We're ready for 3.3.0, which will
hopefully be released soon.

Are the updates for eigen 3.3.0 on GitHub? Looking at the branches, I found one for 3.2.9 but not 3.3.0. 

On Fri, Aug 19, 2016 at 6:07 PM, Bob Carpenter <ca...@alias-i.com> wrote:

It's all done, right?  And isn't it Daniel that figured
out all the traits to make it work?  Not that we're
keeping score.

- Bob

> On Aug 19, 2016, at 11:01 PM, Ben Goodrich <goodri...@gmail.com> wrote:
>
> On Friday, August 19, 2016 at 4:42:55 PM UTC-4, Daniel Lee wrote:
> Bob and Ben have been working through the math library to get it to compile.
>
> Bob has been working. I have been harping.
>
>
> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

> For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

[Charles Margossian]

And to clarify, I mean the updates that make stan/math compatible with eigen 3.3.0. 

[Bob Carpenter]

> On Aug 24, 2016, at 10:16 PM, Charles Margossian <charles.ma...@gmail.com> wrote:
>
> And to clarify, I mean the updates that make stan/math compatible with eigen 3.3.0.

Not sure --- Daniel?

- Bob

[Charles Margossian]

I saw the branch for eigen 3.2.9 on stan-dev/math, but when I went through it, the MatrixExponential class looked the same as in 3.2.8, with most functions and variables being private members. 

I really thought the desired updates were on 3.2.9, but I must have misread their source code. I checked 3.3-beta2, and the MatrixExponential class seems to have the desired revision. 

It probably makes more sense to construct matrix_exp() for eigen 3.3-beta2, rather than fork the class for 3.2.9, and later revise it. That's why I am looking for Stan code updated for 3.3.0.  

On Wed, Aug 24, 2016 at 4:25 PM, Bob Carpenter <ca...@alias-i.com> wrote:

> On Aug 24, 2016, at 10:16 PM, Charles Margossian <charles.margossian93@gmail.com> wrote:
>
> And to clarify, I mean the updates that make stan/math compatible with eigen 3.3.0.

Not sure --- Daniel?

- Bob

[Daniel Lee]

I don't know how close we are to 3.3.0.

Daniel

[Charles Margossian]

I forked the code from 3.3.0 in stan::math, and made some edits to make it compatible with 3.2.9. Prototype code and unit tests are on the branch feature/issue-32-matrix-exponential. Updating when Stan switches to eigen 3.3.0 shouldn't be too hard.

And sorry about the confusion regarding the different versions of eigen... I will be more careful in the future.

On Friday, August 26, 2016 at 11:09:21 AM UTC-4, Daniel Lee wrote:
> I don't know how close we are to 3.3.0.
>
>
>
>
>
>
> Daniel
>
>
>
>
> On Fri, Aug 26, 2016 at 9:45 AM, Charles Margossian <charles.ma...@gmail.com> wrote:
>
> I saw the branch for eigen 3.2.9 on stan-dev/math, but when I went through it, the MatrixExponential class looked the same as in 3.2.8, with most functions and variables being private members. 
>
>
> I really thought the desired updates were on 3.2.9, but I must have misread their source code. I checked 3.3-beta2, and the MatrixExponential class seems to have the desired revision. 
>
>
> It probably makes more sense to construct matrix_exp() for eigen 3.3-beta2, rather than fork the class for 3.2.9, and later revise it. That's why I am looking for Stan code updated for 3.3.0.  
>
>
>
>
> On Wed, Aug 24, 2016 at 4:25 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>
>
> > On Aug 24, 2016, at 10:16 PM, Charles Margossian <charles.ma...@gmail.com> wrote:
>
> >
>
> > And to clarify, I mean the updates that make stan/math compatible with eigen 3.3.0.
>
>
>
> Not sure --- Daniel?
>
>
>
> - Bob
>
>
>
> --
>
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
>

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

>
> For more options, visit https://groups.google.com/d/optout.
>
>
>
>
>
>
>
> --
>
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
>

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.