Repository navigation
Generalize matrix function signatures #1470
Description
Activity
+1 for this. We also need to do this to allow support for sparse matrices. @SteveBronder has been working on that along with perfect forwarding so we can reduce copying and exploit move semantics wherever possible.
We'd need to start coding a bit more defensively after implementing this, otherwise there's the risk of accidentally evaluating the expression multiple times. Take
log_sum_expfor example:template <typename T> auto log_sum_exp(T& x) { if (x.size() == 0) { return -std::numeric_limits<double>::infinity(); } const double max = x.maxCoeff(); if (!std::isfinite(max)) { return max; } return max + std::log((x.array() - max).exp().sum()); }If
log_sum_expis passed an expression, something likelog_sum_exp(m1.array().exp().matrix() * m2), then that expression will get evaluated three times:x.size(),x.maxCoeff(), and the final return statement. So we'd need to start evaluating the expression at the beginning of the function, and probably also updating the documentation so that other/new developers know whyRight. I plan to do that in each function as I generalize its signature. At that time I can also make sure an expression is not used twice. Although I don't think just taking size of an expression evaluates it.
I also agree with the need to document this. Any idea where the doc should be put?
- Absolutely. Is the pattern to do something like: Eigen::Matrix<T::Scalar, -1, 1> x_eval = x.eval() and then use that in the rest of the function? Evaluating `.size()` shouldn't need any evaluation. The expression template will already be evaluated and should at least know its size. `.maxCoeff()` is another story. So we'll have to dig in and see what's costly to evaluate and what's not. It would be a service to devs to update the developer-facing doc. I don't think this calls for any API-level doc changes other than of the acceptable input types for x here in log_sum_exp.…On Nov 29, 2019, at 6:05 PM, Andrew Johnson ***@***.***> wrote: We'd need to start coding a bit more defensively after implementing this, otherwise there's the risk of accidentally evaluating the expression multiple times. Take log_sum_exp for example: template <typename T> auto log_sum_exp(T& x) { if (x.size() == 0) { return -std::numeric_limits<double>::infinity(); } const double max = x.maxCoeff(); if (!std::isfinite(max)) { return max; } return max + std::log((x.array() - max).exp().sum()); } If log_sum_exp is passed an expression, something like log_sum_exp(m1.array().exp().matrix() * m2), then that expression will get evaluated three times: x.size(), x.maxCoeff(), and the final return statement. So we'd need to start evaluating the expression at the beginning of the function, and probably also updating the documentation so that other/new developers know why — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.
- :-) I didn't see this before responding. We're on the same page w.r.t. size(). I think by "evaluate" we mean evaluate operator()(i, j) for every (i, j) pair. It always gets constructed.…On Nov 30, 2019, at 4:00 AM, Tadej Ciglarič ***@***.***> wrote: Right. I plan to do that in each function as I generalize its signature. At that time I can also make sure an expression is not used twice. Although I don't think just taking size of an expression evaluates it. I also agree with the need to document this. Any idea where the doc should be put? — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.
I have a little info on adding new docs to the stan site here
https://github.com/stan-dev/math/blob/develop/doxygen/howto_add_page.md
I wonder if for some functions we use so much info about the matrices it would be better to leave their template arguments as
Matrix<T, R, C>just to force evaluation.If a function needs to evaluate an argument and the template params for that arguments are known (not templated) we can just leave that argument as is -
const Matrix<T, R, C>&. However that does not work if we need template deduction at the same time. In that case I intend to use what Bob suggested.Reacted by Steve BronderIf a function needs to evaluate an argument and the template params for that arguments are known (not templated) we can just leave that argument as is - const Matrix<T, R, C>&
You could also use the
Ref<>functionality to maintain generality of input types, which (AFAIK) will take expressions as inputs and evaluate them into temporaries as needed. Using an example from Eigen's doc:// read-only const argument: void foo2(const Ref<const VectorXf>& x); MatrixXf A; foo2(A.row()); // Compilation error because A.row() is a 1xN object while foo2 is expecting a Nx1 object foo2(A.row().transpose()); // The row is copied into a contiguous temporary foo2(2*a); // The expression is evaluated into a temporary foo2(A.col().segment(2,4)); // No temporaryAlthough some of the copies might be too expensive to be worth it
Reacted by Tadej Ciglarič and Steve BronderFirst, I thought
Refdidn't work with slicing, segmenting, etc., but it seems to work OK. Second, I thought it was type not size that determined passing aVectorXd, but the following compiles and runs.What are the rules for executing into temporaries? I wouldn't have thought any of those template expressions would cause an intermediate.
void foo2(const Eigen::Ref<const Eigen::VectorXd>& x) { std::cout << x << std::endl; } TEST(foo, bar) { Eigen::MatrixXd A(2, 2); A << 1, 2, 3, 4; foo2(A.row(0).transpose()); Eigen::VectorXd a(2); a << 10, 100; foo2(2 * a); foo2(A.col(0).segment(0, 1)); Eigen::MatrixXd B(2, 1); foo2(B); }This page says:
By default, a Ref can reference any dense vector expression of float having a contiguous memory layout. Likewise, a Ref can reference any column-major dense matrix expression of float whose column's elements are contiguously stored with the possibility to have a constant space in-between each column, i.e. the inner stride must be equal to 1, but the outer stride (or leading dimension) can be greater than the number of rows.
In the const case, if the input expression does not match the above requirement, then it is evaluated into a temporary before being passed to the function.
- Thanks---that really clarifies what's going on. I wish we'd figured this out a few years ago!…On Dec 4, 2019, at 2:24 AM, Tadej Ciglarič ***@***.***> wrote: This page says: By default, a Ref can reference any dense vector expression of float having a contiguous memory layout. Likewise, a Ref can reference any column-major dense matrix expression of float whose column's elements are contiguously stored with the possibility to have a constant space in-between each column, i.e. the inner stride must be equal to 1, but the outer stride (or leading dimension) can be greater than the number of rows. In the const case, if the input expression does not match the above requirement, then it is evaluated into a temporary before being passed to the function. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.
I am running out of ideas how to group remaining functions into PRs. So I will just go in alphabetical order.
Reacted by Steve Bronder and Fabio ZotteleI am running out of ideas how to group remaining functions into PRs. So I will just go in alphabetical order.
Can we make a list of the remaining ones as an issue? I think you put one up somewhere a while ago though I couldn't find it. Would be good so we could divvy up work if anyone else is interested. I started picking apart some of the stuff for
read_corr/cov_Letc.No, I never made such a list. If you have a list of all functions in Stan Math I am happy to check ones that I already completed.
I think Steve is reffering to the list of functions we made for the kernel generator.
14 remaining items
It would not get completely banned.
Could we have our
evalfunction handle this?auto x = eval(...);Like if you're using stan::math and you're assigning to an
autovariable, eval. That's how I understand Eigen lol.I guess that means the difference between us and Eigen then is that we can't have expressions as l-values and Eigen can?
Could we have our eval function handle this?
It can be used like this if you want.
I guess that means the difference between us and Eigen then is that we can't have expressions as l-values and Eigen can?
Nope. We are using Eigen expression anyway.
Well Eigen can do:
auto c = (a + b).array().sqrt(); // somehow eigen's sum survives auto d = c + 1.0;And we cannot do:
auto c = softmax(a + b); // our sum gets destructed auto d = c + 1.0;@t4c1 can we start removing the code like
return ret_type(ret);In functions that use
reverse_pass_callback()?https://github.com/stan-dev/math/blob/develop/stan/math/rev/fun/elt_multiply.hpp#L42
If it's safe to pass around arena matrices now I can also put up a lil PR in stanc3 that tests making stuff in the parameters block into arena matrices
Actually wrt stanc3 here's a lil brutal version that uses
arena_tas the matrix type when the scalar type of the matrix is a varstan-dev/stanc3@master...SteveBronder:feature/arena-var-matrix
We could take off the
ret_type()yada from the stuff the uses reverse pass callback rn and get an eyeball on how much we save there. Though imo the bigger thing is going to be doing some form ofarena_t<>with the datacan we start removing the code like
return ret_type(ret);I'm still a bit nervous about this unless something has happened and we're talking about
var<mat>and notmat<var>.Check this code out:
Eigen::Matrix<stan::math::var, -1, -1> myfunc1() { Eigen::MatrixXd res_val(1, 1); res_val << 1.0; stan::arena_t<Eigen::Matrix<stan::math::var, -1, -1>> res = res_val; stan::math::reverse_pass_callback([res]() mutable { std::cout << "myfunc1" << std::endl; std::cout << res.val() << std::endl; }); return res; } auto myfunc2() { Eigen::MatrixXd res_val(1, 1); res_val << 1.0; stan::arena_t<Eigen::Matrix<stan::math::var, -1, -1>> res = res_val; stan::math::reverse_pass_callback([res]() mutable { std::cout << "myfunc2" << std::endl; std::cout << res.val() << std::endl; }); return res; } TEST(AgradRevMatrix, test_arena_output) { auto a = myfunc1(); a(0) = 2.0; auto b = myfunc2(); b(0) = 3.0; stan::math::grad(); }In both cases the function sets the result value to 1, but this is the output:
myfunc2 3 myfunc1 1The
arena_t<Matrix<var, -1, -1>>is playing the same role as thevari **s in the old varis, and if we don't make a copy of that on the output we're in danger of overwriting things.Oof yeah now that I think about this, any
assign()statement could cause this. The only thing I can think of to have it would work be to do the same tricky stuff I'm working on right now forvar<mat>assign where we save the previous slice and re-apply it as a callback in the reverse passWell Eigen can do: And we cannot do:
The difference is that
softmaxneeds to create a local matrix, which is than referenced by returned expression andrsqrtonly needs to reference the given expression.@t4c1 can we start removing the code like
Yeah, Ben is right. We can do it for var. For mat we need the copy. But mat will eventually be gone anyway.
Also I am assuming you are ok with not using holder for arguments.
Reacted by Steve BronderThe difference is that softmax needs to create a local matrix
Yeah and is that not true for all our functions?
I think what we'd write to document this is:
If a Stan expression is going to be used as an lvalue, it shouldn't depend on any rvalue temporaries. Practically this means that an expression saved as an lvalue should only use one function call, because the intermediates allocated in a sequence of function calls will get destructed.
Now that I write that it doesn't sound that clear, but that's what I'm thinking.
I feel like Eigen must have something special going on that we don't cause its expressions merge into bigger and bigger expressions that don't necessarily go away. The reason I'm making the direct comparison to Eigen is it would be cool if we could say "Stan expressions have the same limitations as Eigen expressions", which I don't think is the case cause our functions don't build giant expression types in the same way.
But mat will eventually be gone anyway.
var<mat>design doc needs an update. Assigning to the sub-matrix needs talked about and this seems like a big change of plans too.Also I am assuming you are ok with not using holder for arguments.
Yeah I'm cool with this. It's the same with any oddity. It's fine for the language since we can hide it, and it needs doc'ed for Math so someone doing C++ doesn't get confused and mad.
Reacted by Tadej CiglaričYeah and is that not true for all our functions?
Not for all ,but for many it is true. What I am trying to say is it is not that we are doing things things in different way in our functions - we are doing different things.
it would be cool if we could say "Stan expressions have the same limitations as Eigen expressions"
Actually we are better - we have less limitations! In cases such as softmax that need internal matrices, we will use holder, so that will not be an issue anymore. Eigen's way of doing it would be to simply give up on expressions and eval the result.
I guess softmax was not the best example to start this discussion with, as it will need holder anyway.
Reacted by Ben BalesIn cases such as softmax that need internal matrices, we will use holder, so that will not be an issue anymore
I guess softmax was not the best example to start this discussion with, as it will need holder anyway.Oh okay so we are going to add some holders, just not holders everywhere. Yeah @ me when you put up this pull so I can stare at it a bit.
It will be up in 5 min :)
Reacted by Ben Bales
Description
Matrix functions can be generalized to work with Eigen expressions not just matrices. By propagating matrix expressions around we can avoid creating temporary matrices, improving performance. This could also make adding support for sparse matrices easier.
There was some discoussion about this on discourse. Also an abandoned PR #1309. Contrary to this PR I intend to split the this into smaller chunks:
apply_scalar_unary). Prim implementation of these can in many cases use Eigen implementations that can be faster than std.Example
Instead of :
We can do:
Expected Output
More general and possibly faster code.
Current Version:
v3.0.0