Patch adding matrix/vector scalar subtraction (plus bug fix)
Marcus Brubaker
Also included in the patch is a very small fix to the element-wise division that was recently added. One function was named elt_division when (I assume) it was meant to be elt_divide.
Cheers,
Marcus
Bob Carpenter
one for BEN (ditto).
On 3/7/12 3:26 PM, Marcus Brubaker wrote:
> So I implemented subtraction by a scalar in stan. I've attached a
> patch.
Thanks for the code and heads up about the misnamed
divide operator. Yes, it's very close to the conventions
we've used.
We can add you (Marcus) as a committer to the Stan project
so this doesn't all have to go through e-mails. How do
you feel about doc and unit tests? You can find the tests
under stan/test -- we're using the googletest framework, which
is like JUnit for Java.
?? for EVERYONE:
-----------------
As to this specific patch (and presumably its sibling,
addition by scalar), what does everyone else think about
having this kind of broadcasting of constants up to
appropriately sized vectors or matrices? With v a vector,
m a matrix, and c a scalar, this would be:
vector: (v - c)(i) = v(i) - c
matrix: (m - c)(i,j) = m(i,j) - c
Is this something that's done in MATLAB? In standard
vector notation?
I also have some (long) comments on efficiency and some of
the standards we've tried to follow. Your implementation (for
matrices) is:
inline matrix_d subtract(const matrix_d& m, const double& c) {
matrix_d c_(m.rows(),m.cols());
c_.fill(c);
return m - c_;
}
I'm wondering if it'd be more efficient to do this explicitly:
inline matrix_d subtract(const matrix_d& m, const double c) {
matrix_d mc(m.rows(),m.cols());
for (int i = 0; i < c.size(); ++i)
mc(i) = m(i) - c;
return mc;
}
This unfortunately relies on fixed row-major or column-major
ordering because it uses unary accessors.
But then your last line, "return m - c_;" should be using
a template expression for (m - c_) rather than allocating, so
the return-by-value should copy directly without creating
an intermediate in situations like:
matrix_d mc = m - c;
So they're probably very similar in terms of performance, with
the only difference being the second loop through the data
done by fill(c).
?? for BEN:
-----------
What we'd really like is a no-allocation template expression
for a constant M x N matrix (like your "c_"), but I can't figure
out how to do that in Eigen (CwiseNullaryOp, as returned
by Constant() looks promising). Ben?
You can use size() on matrices and it gives you rows * cols.
This is potentially unsafe if m doesn't have the default
row-major/col-major setting. We've just stuck with the col-major
defaults, so it'll be OK as called from Stan, but not as a general
API function. A full double loop is safer, but much less efficient
in that there's all the index arithmetic
for (int j = 0; j < c.cols(); ++j)
for (int i = 0; i < c.rows(); ++i)
mc(i,j) = m(i,j) - c
We stuck with Eigen's col-major default, so
the loop should at least be col-major.
We relaly want that constant template expression.
We've saved the "x_" style of variables for class member variables
to separate them visually from local variables (we've also been
inconsistent and used "_x" and I used "x__" inside of models to
avoid name clashes with Stan variables).
I haven't profiled, but every piece of advice I've ever
seen says that doubles should be passed by value where
possible. There are some clever ways to template out an arg
so that it's by value for primitives and otherwise by ref --
Alexandrescu discusses this in his C++ book and the approach
is implemented more or less in Boost. But we haven't gone that
far down the optimization path, yet. I'd rather get things
working before getting too carried away on optimization (despite
how I may sound in this e-mail).
We don't want to promote a matrix or vector of scalars to
an auto dif variable (using to_var()) unless we have to,
as the creation of auto-dif variables is expensive.
I plan to go back in and optimize all the vector ops to
avoid many of the to_var() that are there and also partially
evaluate the auto-difs as we've already done for the basic
functions.
The code's essentially the same:
template <typename T1, typename T2>
inline Matrix<var,Dynamic,Dynamic>
subtract(const Matrix<T1,Dynamic,Dynamic>& m,
const T2& c) {
Matrix<var,Dynamic,Dynamic> mc(m.rows(),m.cols());
for (int i = 0; i < c.size(); ++i)
mc(i) = m(i) - c;
return mc;
}
A single templated function could work for all four cases
of double/var for T1 and T2:
template <typename T1, typename T2>
inline Matrix<typename promote_args<T1,T2>::type,Dynamic,Dynamic>
subtract(const Matrix<T1,Dynamic,Dynamic>& m,
const T2& c) {
Matrix<typename promote_args<T1,T2>::type,Dynamic,Dynamic>
mc(m.rows(),m.cols());
for (int i = 0; i < c.size(); ++i)
mc(i) = m(i) - c;
return mc;
}
In fact, if you pull out the two Dynamic template params and replace
them with template params, you could generalize across vectors,
row vectors and matrices, too -- they all support size() and
unary operator() for references.
But as you saw, I've broken these down into separate functions
because all the auto-dif stuff needs a round of unfolding
for optimization (as we've done with the basic functions
you see in stan/agrad/agrad.hpp). And the doc's also a lot
easier than the mess of promotions above.
Also, I've been un-inlining anything that has loops like this.
In general, we only want to inline fairly short programs. Of
course, the compiler probably won't inline the big ones anyway -- it's
only a suggestion.
- Bob
Marcus Brubaker
One question for EVERYONE (marked below as such), and
one for BEN (ditto).Thanks for the code and heads up about the misnamed
On 3/7/12 3:26 PM, Marcus Brubaker wrote:
> So I implemented subtraction by a scalar in stan. I've attached a
> patch.
divide operator. Yes, it's very close to the conventions
we've used.
We can add you (Marcus) as a committer to the Stan project
so this doesn't all have to go through e-mails. How do
you feel about doc and unit tests? You can find the tests
under stan/test -- we're using the googletest framework, which
is like JUnit for Java.
I would be happy to contribute to the doc and unit tests as needed. I'm a bit rusty on my software engineering skills though and have never used either JUnit or googletest so I would ask for you patience. The last serious software engineering (in contrast to research hacking) I've done was probably 10 years ago. I believe JUnit might have just been introduced at the time but maybe not...
?? for EVERYONE:
-----------------
As to this specific patch (and presumably its sibling,
addition by scalar), what does everyone else think about
having this kind of broadcasting of constants up to
appropriately sized vectors or matrices? With v a vector,
m a matrix, and c a scalar, this would be:
vector: (v - c)(i) = v(i) - c
matrix: (m - c)(i,j) = m(i,j) - c
Is this something that's done in MATLAB? In standard
vector notation?
This is MATLAB's convention and to me it is relatively unambiguous and allows for reasonable optimization. The other ways to handle the functionality that I've seen include:
1) require users to define extra "transformed parameters" which manually build vectors and matrices full of the constant
2) provide a ones() function or something similar which produces a matrix (or vector) of appropriate size which can then be multiplied by the desired constant
3) interpret +/- constants as operating on each element of a matrix/vector (what I implemented)
(1) is cumbersome and potentially leads to excess memory usage, although it would be reasonably fast in terms of computation. (2) is the slowest of the options, providing little ability to optimize the code as the intention of the user is hard to discern. That said, functions like ones() and eye() (for identity matrices) are generally handy to have. So I went with (3) since it allows for relatively efficient implementation if done properly (what I did probably isn't).
I also have some (long) comments on efficiency and some of
the standards we've tried to follow. Your implementation (for
matrices) is:
inline matrix_d subtract(const matrix_d& m, const double& c) {
matrix_d c_(m.rows(),m.cols());
c_.fill(c);
return m - c_;
}
I'm wondering if it'd be more efficient to do this explicitly:
inline matrix_d subtract(const matrix_d& m, const double c) {
matrix_d mc(m.rows(),m.cols());
for (int i = 0; i < c.size(); ++i)
mc(i) = m(i) - c;
return mc;
}
This unfortunately relies on fixed row-major or column-major
ordering because it uses unary accessors.
But then your last line, "return m - c_;" should be using
a template expression for (m - c_) rather than allocating, so
the return-by-value should copy directly without creating
an intermediate in situations like:
matrix_d mc = m - c;
So they're probably very similar in terms of performance, with
the only difference being the second loop through the data
done by fill(c).
Assuming Eigen/the compiler is smart, I would think fill() for basic numeric types would reduce to a memset() call which is itself notably faster than looping and repeatedly indexing. I could certainly be wrong.
?? for BEN:
-----------
What we'd really like is a no-allocation template expression
for a constant M x N matrix (like your "c_"), but I can't figure
out how to do that in Eigen (CwiseNullaryOp, as returned
by Constant() looks promising). Ben?
Alternatively, Eigen directly supports adding/subtracting a scalar from each element of an Array type object. Unfortunately, I wasn't able to figure out how to make this work. If I understand Eigen's docs properly you can do it with matrices like so:
return (m.array() - c).matrix();
but I'm not sure about the conversion back to vector types.
We've saved the "x_" style of variables for class member variables
to separate them visually from local variables (we've also been
inconsistent and used "_x" and I used "x__" inside of models to
avoid name clashes with Stan variables).
Ah, sorry, my bad. I'll use that naming convention going forward.
I haven't profiled, but every piece of advice I've ever
seen says that doubles should be passed by value where
possible. There are some clever ways to template out an arg
so that it's by value for primitives and otherwise by ref --
Alexandrescu discusses this in his C++ book and the approach
is implemented more or less in Boost.
In retrospect that's probably true for 64bit machines. On 32bit machines I'd always been told that it's faster to pass by reference since it's half the amount data getting pushed onto the stack. For inlined function it's largely irrelevant, so long as the function doesn't modify the value of the double.
But we haven't gone that
far down the optimization path, yet. I'd rather get things
working before getting too carried away on optimization (despite
how I may sound in this e-mail).
Fair enough. That said I've always found it worth thinking about throughout. Design decisions made now can have performance impacts that could be hard or tedious to fix later.
We don't want to promote a matrix or vector of scalars to
an auto dif variable (using to_var()) unless we have to,
as the creation of auto-dif variables is expensive.
I plan to go back in and optimize all the vector ops to
avoid many of the to_var() that are there and also partially
evaluate the auto-difs as we've already done for the basic
functions.
I'll try to go back and rework what I've done in those terms.
But as you saw, I've broken these down into separate functions
because all the auto-dif stuff needs a round of unfolding
for optimization (as we've done with the basic functions
you see in stan/agrad/agrad.hpp). And the doc's also a lot
easier than the mess of promotions above.
I'm not sure I follow you in terms of the docs. Having dozens of "add" functions compared to one seems more confusing to me or am I missing your meaning?
Also, I've been un-inlining anything that has loops like this.
In general, we only want to inline fairly short programs. Of
course, the compiler probably won't inline the big ones anyway -- it's
only a suggestion.
Yeah, I've generally understood inline to be more about allowing the compiler to inline it if it chooses to, rather than forcing it to do anything. My only argument for keeping those inlined is that ideally most of those loops would eventually go away and be replaced by modular calls to other appropriate functions and at that point someone would need to remember to add the inline designation back in.
Cheers,
Marcus
Marcus Brubaker
I realized that the various vector types are just specialized Matrix templates (a retrospective duh) so I can use the .array()/.matrix() trick to implement the addition/subtraction of elements cleanly and efficiently. I've also templated the functions to avoid code duplication, removed the pass-by-reference doubles and implemented addition.
By the way, is there a reason to differentiate the matrix/vector/row_vector functions instead of templating them?
I've attached a patch for now. If I get commit access I can commit it assuming that it passes muster.
I think the .array()/.matrix() trick could be used to avoid the loops in the elementwise multiplication and division as well as I think those operators are defined there. If you'd like I can go ahead and do that and send a patch.
Cheers,
Marcus
Andrew Gelman
On Mar 7, 2012, at 5:06 PM, Bob Carpenter wrote:
>
> ?? for EVERYONE:
> -----------------
> As to this specific patch (and presumably its sibling,
> addition by scalar), what does everyone else think about
> having this kind of broadcasting of constants up to
> appropriately sized vectors or matrices? With v a vector,
> m a matrix, and c a scalar, this would be:
>
> vector: (v - c)(i) = v(i) - c
>
> matrix: (m - c)(i,j) = m(i,j) - c
>
This looks good to me.
Bob Carpenter
in my mailer (Thunderbird). I've been meaning to upgrade
to something else anyway now that I'm on a Mac.
I have the ">" convention going, so I'll comment inline.
On 3/7/12 6:14 PM, Marcus Brubaker wrote:
> On Wed, Mar 7, 2012 at 5:06 PM, Bob Carpenter <ca...@alias-i.com
> <mailto:ca...@alias-i.com>> wrote:
>
> One question for EVERYONE (marked below as such), and
> one for BEN (ditto).
>
>
> On 3/7/12 3:26 PM, Marcus Brubaker wrote:
> > So I implemented subtraction by a scalar in stan. I've attached a
> > patch.
>
> Thanks for the code and heads up about the misnamed
> divide operator. Yes, it's very close to the conventions
> we've used.
>
> We can add you (Marcus) as a committer to the Stan project
> so this doesn't all have to go through e-mails. How do
> you feel about doc and unit tests? You can find the tests
> under stan/test -- we're using the googletest framework, which
> is like JUnit for Java.
>
>
> I would be happy to contribute to the doc and unit tests as needed. I'm
> a bit rusty on my software engineering skills though and have never used
> either JUnit or googletest so I would ask for you patience. The last
> serious software engineering (in contrast to research hacking) I've done
> was probably 10 years ago. I believe JUnit might have just been
> introduced at the time but maybe not...
JUnit was definitely around ten years ago. But the concept's
older. The framework's pretty simple if you look at the tests
and Daniel has everything instrumented so you can just run tests
by:
% cd $STAN
% make test/math/matrix
The test itself is in $STAN/src/test/math/matrix_test.cpp.
There's a matching test in test/agrad/matrix for the auto-dif
stuff.
> ?? for EVERYONE:
> -----------------
> As to this specific patch (and presumably its sibling,
> addition by scalar), what does everyone else think about
> having this kind of broadcasting of constants up to
> appropriately sized vectors or matrices? With v a vector,
> m a matrix, and c a scalar, this would be:
>
> vector: (v - c)(i) = v(i) - c
>
> matrix: (m - c)(i,j) = m(i,j) - c
>
> Is this something that's done in MATLAB? In standard
> vector notation?
>
>
>
> This is MATLAB's convention and to me it is relatively unambiguous and
> allows for reasonable optimization.
I don't see any problem with ambiguity. And it'll definitely
be the most optimizable approach.
> The other ways to handle the
> functionality that I've seen include:
>
> 1) require users to define extra "transformed parameters" which manually
> build vectors and matrices full of the constant
> 2) provide a ones() function or something similar which produces a
> matrix (or vector) of appropriate size which can then be multiplied by
> the desired constant
> 3) interpret +/- constants as operating on each element of a
> matrix/vector (what I implemented)
>
> (1) is cumbersome and potentially leads to excess memory usage, although
> it would be reasonably fast in terms of computation.
Actually, (1) is independent of (2) and (3). We could construct
ahead the templated constant and re-use it, which would be even
more efficient than building it each time (though probably close, it's
a very simple function).
(3) is definitely what we want, but I we'll get it through
template expressions. At which point we can optimize further
by pre-building the template expression as in (1).
If we have a const(M,N,x) function, we can expose (2) and get
the efficiency of (3). That's the beauty of template expressions.
> (2) is the slowest
> of the options, providing little ability to optimize the code as the
> intention of the user is hard to discern. That said, functions like
> ones() and eye() (for identity matrices) are generally handy to have.
I don't think I could bear the name "eye()". Tell me they don't
use that in MATLAB. Please. People used "liszt" in Prolog because
"list" was reserved, and it drove me bonkers. They use "clazz" in
Java because "class" is reserved.
You couldn't have a ones() function as is, it'd have to be ones(M,N)
with the sizes. This is actually what Eigen does, though I now see
they implement it efficiently storing only the sizes and value.
> So I went with (3) since it allows for relatively efficient
> implementation if done properly (what I did probably isn't).
The best implementation would be to have a very small on-the-stack
template expression, and it turns out that's what Eigen does (see
below).
I finally bit the bullet and tracked down what was
going on in Eigen (I was more curious than concerned about
the case we're looking at).
Your intuitions are right, but for the wrong reason.
fill() is just an alias for setConstant, which does this
(from CwiseNullaryOp.h in Eigen/Core):
template<typename Derived>
EIGEN_STRONG_INLINE Derived&
PlainObjectBase<Derived>::setConstant(Index rows, Index cols, const
Scalar& value)
{
resize(rows, cols);
return setConstant(value);
}
which calls this:
template<typename Derived>
EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::setConstant(const
Scalar& value)
{
return derived() = Constant(rows(), cols(), value);
}
The derived() thing is from their use of the so-called curiously
recurring template pattern. This then calls this method:
template<typename Derived>
EIGEN_STRONG_INLINE const typename DenseBase<Derived>::ConstantReturnType
DenseBase<Derived>::Constant(Index rows, Index cols, const Scalar& value)
{
return DenseBase<Derived>::NullaryExpr(rows, cols,
internal::scalar_constant_op<Scalar>(value));
}
which then calls:
template<typename Derived>
template<typename CustomNullaryOp>
EIGEN_STRONG_INLINE const CwiseNullaryOp<CustomNullaryOp, Derived>
DenseBase<Derived>::NullaryExpr(Index rows, Index cols, const
CustomNullaryOp& func)
{
return CwiseNullaryOp<CustomNullaryOp, Derived>(rows, cols, func);
}
and finally, we construct the CwiseNullaryOp, which only stores the
number of rows, number of columns, and the functor (which is just
a single value promoted in the second to last call above).
> ?? for BEN:
> -----------
> What we'd really like is a no-allocation template expression
> for a constant M x N matrix (like your "c_"), but I can't figure
> out how to do that in Eigen (CwiseNullaryOp, as returned
> by Constant() looks promising). Ben?
>
>
> Alternatively, Eigen directly supports adding/subtracting a scalar from
> each element of an Array type object. Unfortunately, I wasn't able to
> figure out how to make this work. If I understand Eigen's docs properly
> you can do it with matrices like so:
>
> return (m.array() - c).matrix();
>
> but I'm not sure about the conversion back to vector types.
Ben liked the array() thingy too. I can look into playing
nicely with vectors.
> We've saved the "x_" style of variables for class member variables
> to separate them visually from local variables (we've also been
> inconsistent and used "_x" and I used "x__" inside of models to
> avoid name clashes with Stan variables).
>
>
> Ah, sorry, my bad. I'll use that naming convention going forward.
I think that's more our inconsistency.
> I haven't profiled, but every piece of advice I've ever
> seen says that doubles should be passed by value where
> possible. There are some clever ways to template out an arg
> so that it's by value for primitives and otherwise by ref --
> Alexandrescu discusses this in his C++ book and the approach
> is implemented more or less in Boost.
>
>
> In retrospect that's probably true for 64bit machines. On 32bit
> machines I'd always been told that it's faster to pass by reference
> since it's half the amount data getting pushed onto the stack. For
> inlined function it's largely irrelevant, so long as the function
> doesn't modify the value of the double.
That's the general advice for both 32-bit and 64-bit situations.
It will probably depend on how the variable gets used in the function,
but the problem with refs is that you need to chase pointers, then
bring the 64 bits into registers. This is going to be
way more expensive than copying an extra 32 bits in terms
of memory non-locality, which is the main killer given how much
faster CPUs are than memory these days.
Even some structs, like our var, could probably be as efficiently,
if not more efficiently, passed by value. They're just pointers
to implementation, so they're only 8 bytes.
Again, this is something we can play with later.
> But we haven't gone that
> far down the optimization path, yet. I'd rather get things
> working before getting too carried away on optimization (despite
> how I may sound in this e-mail).
>
>
> Fair enough. That said I've always found it worth thinking about
> throughout. Design decisions made now can have performance impacts that
> could be hard or tedious to fix later.
I agree. The issue for me is that they need to be designed
so that they can be optimized later. In that sense, we can
add (v - c), and worry about how to implement it more efficiently
later without changing the modeling language. That's all I'm
really worried about now.
> We don't want to promote a matrix or vector of scalars to
> an auto dif variable (using to_var()) unless we have to,
> as the creation of auto-dif variables is expensive.
> I plan to go back in and optimize all the vector ops to
> avoid many of the to_var() that are there and also partially
> evaluate the auto-difs as we've already done for the basic
> functions.
>
>
> I'll try to go back and rework what I've done in those terms.
>
> But as you saw, I've broken these down into separate functions
> because all the auto-dif stuff needs a round of unfolding
> for optimization (as we've done with the basic functions
> you see in stan/agrad/agrad.hpp). And the doc's also a lot
> easier than the mess of promotions above.
>
>
> I'm not sure I follow you in terms of the docs. Having dozens of "add"
> functions compared to one seems more confusing to me or am I missing
> your meaning?
It's that each one has simpler doc that doesn't involve
calculated return types.
But I doubt anyone's going to use the doxygen doc anyway.
It's almost useless for Boost and Eigen in my experience,
given their insane levels of complexity.
> Also, I've been un-inlining anything that has loops like this.
> In general, we only want to inline fairly short programs. Of
> course, the compiler probably won't inline the big ones anyway -- it's
> only a suggestion.
>
>
> Yeah, I've generally understood inline to be more about allowing the
> compiler to inline it if it chooses to, rather than forcing it to do
> anything. My only argument for keeping those inlined is that ideally
> most of those loops would eventually go away and be replaced by modular
> calls to other appropriate functions and at that point someone would
> need to remember to add the inline designation back in.
Either way, we'll profile when we want to do the next
speed improvements, and that'll tell us where the
hot spots are in practice.
- Bob
Marcus Brubaker
Oddly your replies don't show up properly formatted
in my mailer (Thunderbird). I've been meaning to upgrade
to something else anyway now that I'm on a Mac.
Strange. I use thunderbird as well for my other accounts. I'm not sure why they'd show up differently but things look fine from my end with Gmail.
JUnit was definitely around ten years ago. But the concept's
older. The framework's pretty simple if you look at the tests
and Daniel has everything instrumented so you can just run tests
by:
% cd $STAN
% make test/math/matrix
The test itself is in $STAN/src/test/math/matrix_test.cpp.
There's a matching test in test/agrad/matrix for the auto-dif
stuff.
I'll take a look and try to add a few tests where appropriate.
The other ways to handle the
functionality that I've seen include:
1) require users to define extra "transformed parameters" which manually
build vectors and matrices full of the constant
2) provide a ones() function or something similar which produces a
matrix (or vector) of appropriate size which can then be multiplied by
the desired constant
3) interpret +/- constants as operating on each element of a
matrix/vector (what I implemented)
Actually, (1) is independent of (2) and (3). We could construct
ahead the templated constant and re-use it, which would be even
more efficient than building it each time (though probably close, it's
a very simple function).
True, although unless the modeling language (or compiler) recognizes a matrix of constants it would be hard to prevent it from falling back to performing the operation as on two full matrices/vectors.
(3) is definitely what we want, but I we'll get it through
template expressions. At which point we can optimize further
by pre-building the template expression as in (1).
If we have a const(M,N,x) function, we can expose (2) and get
the efficiency of (3). That's the beauty of template expressions.
Agreed.
I don't think I could bear the name "eye()". Tell me they don't
(2) is the slowest
of the options, providing little ability to optimize the code as the
intention of the user is hard to discern. That said, functions like
ones() and eye() (for identity matrices) are generally handy to have.
use that in MATLAB. Please. People used "liszt" in Prolog because
"list" was reserved, and it drove me bonkers. They use "clazz" in
Java because "class" is reserved.
Sadly, yes, there is an eye function in Matlab. I've been using it too long now, but I remember finding it idiotic when I first started with Matlab. It's gotten so bad that in Matlab dominated fields some papers will use it as mathematical notation. id() would be a better option.
You couldn't have a ones() function as is, it'd have to be ones(M,N)
with the sizes. This is actually what Eigen does, though I now see
they implement it efficiently storing only the sizes and value.
Sorry, I was a little sloppy and unclear. I assumed there to be size parameters passed, I was just being lazy.
The best implementation would be to have a very small on-the-stack
So I went with (3) since it allows for relatively efficient
implementation if done properly (what I did probably isn't).
template expression, and it turns out that's what Eigen does (see
below).
Agreed. I think the .array()/.matrix() approach does this quite well and with a bit less obfuscation although I haven't dug into what, exactly, Eigen is doing with it. According to it's docs .array() and .matrix() should be effectively free once the compiler does its optimization.
I finally bit the bullet and tracked down what was
Assuming Eigen/the compiler is smart, I would think fill() for basic
numeric types would reduce to a memset() call which is itself notably
faster than looping and repeatedly indexing. I could certainly be wrong.
going on in Eigen (I was more curious than concerned about
the case we're looking at).
Your intuitions are right, but for the wrong reason.
You are a braver man than I chasing that down... But it's very comforting to know that Eigen is doing something rather sensible in that case.
Ben liked the array() thingy too. I can look into playing
?? for BEN:
-----------
What we'd really like is a no-allocation template expression
for a constant M x N matrix (like your "c_"), but I can't figure
out how to do that in Eigen (CwiseNullaryOp, as returned
by Constant() looks promising). Ben?
Alternatively, Eigen directly supports adding/subtracting a scalar from
each element of an Array type object. Unfortunately, I wasn't able to
figure out how to make this work. If I understand Eigen's docs properly
you can do it with matrices like so:
return (m.array() - c).matrix();
but I'm not sure about the conversion back to vector types.
nicely with vectors.
It turns out that it works just fine with vectors, I stopped digging too soon. See my most recent patch. I think this actually opens up a number of nice implementation approaches for these kinds of operators because the Array class has a bunch of operations that Matrix doesn't.
That's the general advice for both 32-bit and 64-bit situations.
I haven't profiled, but every piece of advice I've ever
seen says that doubles should be passed by value where
possible. There are some clever ways to template out an arg
so that it's by value for primitives and otherwise by ref --
Alexandrescu discusses this in his C++ book and the approach
is implemented more or less in Boost.
In retrospect that's probably true for 64bit machines. On 32bit
machines I'd always been told that it's faster to pass by reference
since it's half the amount data getting pushed onto the stack. For
inlined function it's largely irrelevant, so long as the function
doesn't modify the value of the double.
It will probably depend on how the variable gets used in the function,
but the problem with refs is that you need to chase pointers, then
bring the 64 bits into registers. This is going to be
way more expensive than copying an extra 32 bits in terms
of memory non-locality, which is the main killer given how much
faster CPUs are than memory these days.
Very true. I've dropped the pass-by-reference in my latest patch.
It's that each one has simpler doc that doesn't involve
I'm not sure I follow you in terms of the docs. Having dozens of "add"
functions compared to one seems more confusing to me or am I missing
your meaning?
calculated return types.
I don't really see the calculated return types as an issue. Personally, I find multiple versions of the same function distracting as opposed to a single templated function. The added advantage of the single templated function is improved testing coverage as there are fewer redundant code paths to check.
But I doubt anyone's going to use the doxygen doc anyway.
It's almost useless for Boost and Eigen in my experience,
given their insane levels of complexity.
Sad but true. Just today I was trying figure out something with Eigen and their doxygen docs were next to useless even though they are fairly complete for precisely this reason.
Cheers,
Marcus
Andrew Gelman
I don't think I could bear the name "eye()". Tell me they don't
use that in MATLAB. Please. People used "liszt" in Prolog because
"list" was reserved, and it drove me bonkers. They use "clazz" in
Java because "class" is reserved.
Sadly, yes, there is an eye function in Matlab. I've been using it too long now, but I remember finding it idiotic when I first started with Matlab. It's gotten so bad that in Matlab dominated fields some papers will use it as mathematical notation. id() would be a better option.
Marcus Brubaker
Ben Goodrich
I realized that the various vector types are just specialized Matrix templates (a retrospective duh) so I can use the .array()/.matrix() trick to implement the addition/subtraction of elements cleanly and efficiently. I've also templated the functions to avoid code duplication, removed the pass-by-reference doubles and implemented addition.
Just sayin' that the .array() / operate / .matrix() route is also the way I would have done it.
P.S. @Bob, note there is an git apply function for patches.