Using multiple algorithms for one function

[Charles Margossian]

Related to this issue (https://github.com/stan-dev/math/issues/32#issuecomment-244234095), but seems general enough to be on the dev list.

I have multiple algorithms for computing the matrix exponential. A general algorithm, and then some specialized ones that target specific kinds of matrices.

Should the algorithms be saved in different files? Right now, I have two files: one with the pade approximation from Eigen 3.0, and another one with code for symmetric, nilpotent, and 2x2 matrices. Should they each be stored in separate files? Should each algorithms have its own unit tests, or is testing the final expm() function sufficient?

I'm also unclear on when to write a template function, and when to overload it for fvar<T> and var. Doing the latter makes sense to me, if we handwrite the gradient (in the fvar version), else would it be redundant?

[Bob Carpenter]

> On Sep 6, 2016, at 1:33 PM, Charles Margossian <charles.ma...@gmail.com> wrote:
>
> Related to this issue (https://github.com/stan-dev/math/issues/32#issuecomment-244234095), but seems general enough to be on the dev list.
>
> I have multiple algorithms for computing the matrix exponential. A general algorithm, and then some specialized ones that target specific kinds of matrices.

What are the function names going to be? For each function,
you want a separate file for double, rev, and fvar<T>:

* double-only versions should go in prim/mat/fun
* var versions in rev/mat/fun
* fvar<T> template functions in fwd/mat/fun

Basically, you should follow the pattern for things like
the Cholesky decomposition (I think that has custom
gradients).

> Should the algorithms be saved in different files? Right now, I have two files: one with the pade approximation from Eigen 3.0, and another one with code for symmetric, nilpotent, and 2x2 matrices.

The 2 x 2 case should be rolled into the general function.

Symmetric matrix input goes in a different function, I think,
though Ben seemed to think you should just overload a single
function and choose the algorithm.

> Should they each be stored in separate files? Should each algorithms have its own unit tests, or is testing the final expm() function sufficient?

Yes, every instanatiation of every function needs to
be tested.

You'll probably at least have:

expm
expm_nilpotent

and maybe

expm_symmetric

if that doesn't get rolled into expm.

>
> I'm also unclear on when to write a template function, and when to overload it for fvar<T> and var.

fvar<T> requires a template function. Others don't.

> Doing the latter makes sense to me, if we handwrite the gradient (in the fvar version), else would it be redundant?

You want to make the var version efficient. Again, you'll
want to follow the other examples of matrix derivatives.

Rob may be the best person to help with this at this point.

- 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.

[Charles Margossian]

What are the function names going to be?

I plan to create a single function: expm(). The function checks for certain properties of the input matrix and then picks an algorithm. 

I've written expm() under prim/mat/fun, rev/mat/fun, and fwd/mat/fun, which then calls other functions, such as matrix_exp_compute_2x2. Currently, I didn't handwrite the gradients. 

matrix_exp_compute_2x2 is stored in prim/mat/fun and has the following declaration:

template <typename T>

        void matrix_exp_compute_2x2(const Matrix<T, Dynamic, Dynamic>& A,

                               Matrix<T, Dynamic, Dynamic>& B)

 

Because it is a template type, it properly handles T = var, T = fvar<T'>, and T = double. I don't see the point of writing the function three times for var, fvar<T>, and double. Same rule of thumb applies for matrix_exp_compute_sym and matrix_exp_compute_nilpotent. 

You want to make the var version efficient.  Again, you'll
want to follow the other examples of matrix derivatives.

 

Ok, this is where I see the point of overloading the function for var and fvar<T>. Clearly, I'm not making the var version efficient yet. This is where I need to dig. 

Yes, every instantiation of every function needs to
be tested.

I have unit tests that call expm() and test it on different matrices so that every algorithm gets tested. I don't think having a unit test that explicitly calls matrix_exp_compute_2x2 is needed, since users will call expm().   

On Tue, Sep 6, 2016 at 7:40 PM, Bob Carpenter <ca...@alias-i.com> wrote:

> 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.

[Bob Carpenter]

> On Sep 7, 2016, at 10:48 AM, Charles Margossian <charles.ma...@gmail.com> wrote:
>
> What are the function names going to be?
>
> I plan to create a single function: expm(). The function checks for certain properties of the input matrix and then picks an algorithm.

How easy it to check nilpotency? I'm OK with symmetry test even
if it'll slow down non-symmetric calls a bit, but it's tricky because
matrices are sometimes almost symmetric, but not technically symmetric.
If you use ==, you'll almost never call it, and if you don't, you have
to decide on an absolute or relative symmetry tolerance.
An absolute check would depend on scale, but the check could be made
relative at quite a bit more cost (two divisions and an extra sum per entry).

If you havne't specialized, relative symmetry tests are going to be expensive
as they will add to the autodiff stack.

> I've written expm() under prim/mat/fun, rev/mat/fun, and fwd/mat/fun, which then calls other functions, such as matrix_exp_compute_2x2. Currently, I didn't handwrite the gradients.

Then how did you implement in rev and fwd? Those reqire
gradients to be coded. Otherwise, you only need the templated
version in prim.

>
> matrix_exp_compute_2x2 is stored in prim/mat/fun and has the following declaration:
>
> template <typename T>
>
> void matrix_exp_compute_2x2(const Matrix<T, Dynamic, Dynamic>& A,
>
> Matrix<T, Dynamic, Dynamic>& B)

Why are you including a mutable matrix B as an argument rather than
just returning a matrix?

> Because it is a template type, it properly handles T = var, T = fvar<T'>, and T = double. I don't see the point of writing the function three times for var, fvar<T>, and double. Same rule of thumb applies for matrix_exp_compute_sym and matrix_exp_compute_nilpotent.

The point is that then you can code derivatives directly, which
is much much more efficient than autodiffing in both memory and
speed.

> You want to make the var version efficient. Again, you'll
> want to follow the other examples of matrix derivatives.
>
> Ok, this is where I see the point of overloading the function for var and fvar<T>. Clearly, I'm not making the var version efficient yet. This is where I need to dig.

Ben pointed to the Giles paper's pointers about how to calculate
the derivatives explicitly. The point is not just efficiency, but
that autodiffing the approximation algorithm for expm() isn't necessarily
going to give you good approximations to the derivatives. And it
won't be efficient.

> Yes, every instantiation of every function needs to
> be tested.
>
> I have unit tests that call expm() and test it on different matrices so that every algorithm gets tested. I don't think having a unit test that explicitly calls matrix_exp_compute_2x2 is needed, since users will call expm().

The reason to test it separately is that otherwise your tests
for expm() need to know the branching logic internally to perform
adequate testing. Daniel will probably have an opinion on this, too.

- Bob

>
>
> On Tue, Sep 6, 2016 at 7:40 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>

> > 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.

> 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.

[Charles Margossian]

How easy it to check nilpotency?

Not sure how to best do this. The check I wrote involves an embedded for loop and uses ==. I might hold off doing the nilpotent specialization for now. I know matlab handles it, but I'm not aware of any concrete applications (at least not in pharmacometrics). It's not hard per se -- I have a working prototype --, but computing a Taylor series with matrices probably requires more careful coding, at least to be optimized. 

Then how did you implement in rev and fwd?  Those reqire
gradients to be coded.  Otherwise, you only need the templated
version in prim.

Ok, that was my main point of confusion and it now makes sense! So yes, I'll code the gradients, following Giles' paper, in the fvar and var versions. I also didn't realize how "risky" relying on auto-diff approximation could be. 

Why are you including a mutable matrix B as an argument rather than
just returning a matrix? 

I am being consistent with matrix_exp_compute() in Eigen/unsupported, but I don't think there is any advantage in doing so. I'll switch to a simpler signature. In Eigen it makes sense because of the way they store the result of A.exp() -- the solution B is a class member. 

The reason to test it separately is that otherwise your tests
for expm() need to know the branching logic internally to perform
adequate testing.  Daniel will probably have an opinion on this, too.

I see. I propose creating multiple tests in one unit test file then, to indicate these algorithms all ultimately serve the same purpose (computing a matrix exponential). I'm happy to hear Daniel's opinion on the matter.  

- Charles 

On Wed, Sep 7, 2016 at 11:00 AM, Bob Carpenter <ca...@alias-i.com> wrote:

> > 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.

> 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.

> 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.

[Avraham Adler]

On Wednesday, September 7, 2016 at 11:55:09 AM UTC-4, Charles Margossian wrote:
> How easy it to check nilpotency?
>
>
> Not sure how to best do this. The check I wrote involves an embedded for loop and uses ==. I might hold off doing the nilpotent specialization for now. I know matlab handles it, but I'm not aware of any concrete applications (at least not in pharmacometrics). It's not hard per se -- I have a working prototype --, but computing a Taylor series with matrices probably requires more careful coding, at least to be optimized. 

Excuse the ignoramus intrusion, but per https://en.wikipedia.org/wiki/Nilpotent_matrix#Characterization, iff the only eigenvalue of a square matrix is 0 it is nilpotent. Would it be computationally more expensive to calculate the eigenvectors?

Avi

[Bob Carpenter]

> On Sep 7, 2016, at 11:55 AM, Charles Margossian <charles.ma...@gmail.com> wrote:
>
> How easy it to check nilpotency?
>
> Not sure how to best do this. The check I wrote involves an embedded for loop and uses ==.

Loops and == don't do any allocation on the autodiff stack, but
they're not so efficient if the for loop is big, which is why
you'd want to put them in a separate function. Presumably the users
know when they're passing in these kinds of results.

My general advice is to get a version working first, even if
it's just one templated implementation in prim/mat and then
tests in prim/mat, prim/rev, and prim/fwd. Then you can build
out making things more efficient. And definitely don't prioritize
things we don't have a use for immediately.

> I might hold off doing the nilpotent specialization for now. I know matlab handles it, but I'm not aware of any concrete applications (at least not in pharmacometrics). It's not hard per se -- I have a working prototype --, but computing a Taylor series with matrices probably requires more careful coding, at least to be optimized.
>
> Then how did you implement in rev and fwd? Those reqire
> gradients to be coded. Otherwise, you only need the templated
> version in prim.
>
> Ok, that was my main point of confusion and it now makes sense! So yes, I'll code the gradients, following Giles' paper, in the fvar and var versions. I also didn't realize how "risky" relying on auto-diff approximation could be.

The specific risks are numerical accuracy and efficiency
both of which can be tested, but we can only test numerical
accuracy if we have a way of creating a more accurate result.

You can also test with finite differences, though we know those
won't be very accurate.

> Why are you including a mutable matrix B as an argument rather than
> just returning a matrix?
>
> I am being consistent with matrix_exp_compute() in Eigen/unsupported, but I don't think there is any advantage in doing so. I'll switch to a simpler signature. In Eigen it makes sense because of the way they store the result of A.exp() -- the solution B is a class member.

It would also make sense if you'd already allocated
a result and then were calling this function to fill it in.
So that was a real question, not a rhetorical one. I
couldn't tell from just the fragment.

> I see. I propose creating multiple tests in one unit test file then, to indicate these algorithms all ultimately serve the same purpose (computing a matrix exponential). I'm happy to hear Daniel's opinion on the matter.

Each function should go in its own file, and primitive, reverse,
and forward versions all go in their own directories.

Tests for the functions should follow the same file structure.

- Bob

[Bob Carpenter]

We definitely don't want to perform a general eigendecomposition
to test whether we can apply a more efficient matrix exponential
algorithm. Is there a fast way to just find the largest eigenvalue?
If not, you'd want to break the nilpotent case out. I have no
idea what it's used for and if you don't have an application in
mind, I'd just recommend skipping it for now.

- 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]

My general advice is to get a version working first, even if
it's just one templated implementation in prim/mat and then
tests in prim/mat, prim/rev, and prim/fwd.  

Should I make a pull request once I have the simple version in place? 

> 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 then create a new issue for optimizing)?

[Daniel Lee]

Yes! Definitely work in stages. Having a working version with tests makes it very easy to replace with better code later.

On Thu, Sep 8, 2016 at 8:39 AM, Charles Margossian <charles.ma...@gmail.com> wrote:

(and then create a new issue for optimizing)?

[Krzysztof Sakrejda]

Yes to both! This whole process works better with small well tested pull requests. . K

On Thu, Sep 8, 2016, 8:39 AM Charles Margossian <charles.ma...@gmail.com> wrote:

(and then create a new issue for optimizing)?

On Thu, Sep 8, 2016 at 8:38 AM, Charles Margossian <charles.ma...@gmail.com> wrote:

My general advice is to get a version working first, even if
it's just one templated implementation in prim/mat and then
tests in prim/mat, prim/rev, and prim/fwd.  

Should I make a pull request once I have the simple version in place? 

On Wed, Sep 7, 2016 at 7:55 PM, Bob Carpenter <ca...@alias-i.com> wrote:

We definitely don't want to perform a general eigendecomposition
to test whether we can apply a more efficient matrix exponential
algorithm.  Is there a fast way to just find the largest eigenvalue?
If not, you'd want to break the nilpotent case out.  I have no
idea what it's used for and if you don't have an application in
mind, I'd just recommend skipping it for now.

- Bob

> On Sep 7, 2016, at 6:04 PM, Avraham Adler <avraha...@gmail.com> wrote:
>
> On Wednesday, September 7, 2016 at 11:55:09 AM UTC-4, Charles Margossian wrote:
>> How easy it to check nilpotency?
>>
>>
>> Not sure how to best do this. The check I wrote involves an embedded for loop and uses ==. I might hold off doing the nilpotent specialization for now. I know matlab handles it, but I'm not aware of any concrete applications (at least not in pharmacometrics). It's not hard per se -- I have a working prototype --, but computing a Taylor series with matrices probably requires more careful coding, at least to be optimized.
>
>
> Excuse the ignoramus intrusion, but per https://en.wikipedia.org/wiki/Nilpotent_matrix#Characterization, iff the only eigenvalue of a square matrix is 0 it is nilpotent. Would it be computationally more expensive to calculate the eigenvectors?
>
> Avi
>
> --
> 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.

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.

[Daniel Lee]

Regarding testing an internal function: yes, please test and document it. Please.

I can't stress enough the number of times we hit boundary conditions where the internal functions aren't so well behaved with the inputs that its given and there aren't any tests. It takes a lot more time and energy then to even figure out how a bit of code that was written some time ago can be instantiated, let alone figure out what it's supposed to do.

So yes, test the internal function. By test, I think people are thinking it's a large burden, but what I'm looking for are essentially these things:

- existence of tests. Getting over the hurdle from 0 to 1 test is huge. Coming back to the code and having a test there waiting for you to extend is so much easier than trying from scratch.

- test that it works well within the range of input values that you're expecting. That could be one test. It could be many. But if it's going out to the public (even through another function), we should at least check once that it behaves well.

- test that it fails the way you intend it to fail. In some cases, especially if the function is well doced, input that violates the preconditions will be undefined behavior. Show that in a test. Or in other cases, input that is outside some range should throw or do something else. Check that.

If you're writing some code, you should do at least those things. Then there's branch coverage, code coverage, etc. which I don't really expect unless it's clear that there's completely different behavior.

So yes, please test internal functions. Especially if they're approximations that only hold for a particular range of values.

[Krzysztof Sakrejda]

It also makes sense to do this testing at the lowest level possible so if you have a higher function that choose among approximations you'll only be testing that it's branching correctly rather than testing the complete set of algorithms.

> > 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.

> 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.

> 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.

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.

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.

[Daniel Lee]

+1

[Bob Carpenter]

Just a programming point here that slightly contradicts what
Daniel was saying

Internal functions should *not* check their input conditions.
So you don't want the symmetric handler to check for symmetry
or the 2x2 handler to check sizes. That's the job of the calling
function and should be tested from there. It helps ensure the
error reports are correct (they come from the function the user
called, not some internal function they've never heard of)
and removes redundancy in tests.

You can test that internal errors arise from the wrong inputs,
but I wouldn't bother.

There are two drawbacks to testing internal functions:

* when the internal functions change without changing the external
API, then the tests need to change even though it doesn't
change external public behavior

* requiring tests for internal functions discourages developers
from breaking natural functional units into tests

This would all be mitigated somewhat if the testing framework
were faster. Right now, with a lot of functions not just
including what they use, a change in a low-level functions requires
a 45 minute unit test run on my end, or I often find failures
on Jenkins, and have to go back and do everything again.

- Bob

> > > 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.

> > 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.

> > 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.

> 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.

> 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]

-1, since we're voting :-)

The top-level function *must* test that it's branching correctly
and that the branches return the right results. It also *must* test
that it returns appropriate error conditions everywhere.

- Bob

[Krzysztof Sakrejda]

On Thu, Sep 8, 2016 at 10:55 AM Bob Carpenter <ca...@alias-i.com> wrote:

-1, since we're voting :-)

The top-level function *must* test [snip] that the branches return the right results.  

If the top level function tests that it's branching correctly and the lower-level functions test that they return the correct results then having the top level function test for correct results only duplicates tests. Maybe ideally we'd love to have the top-level functions do that but if the top-level function test fails you have no idea what level the calculation failed at.  We were just trying to hash this out w.r.t. the tests for the Stan refactor and I thought the conclusion was that the function should test it's own code and assume the lower levels are tested separately... 

[Daniel Lee]

+1 to all of Bob's comments. Especially the part about internal functions and inputs. They shouldn't check so you don't need to test inputs. But doc them! Please.

>> > > 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.

>> > 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.

>> > 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.

>> 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.

>> 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.

>> 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.

>> 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.

> 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.

[Krzysztof Sakrejda]

On Thu, Sep 8, 2016 at 10:54 AM Bob Carpenter <ca...@alias-i.com> wrote:

Just a programming point here that slightly contradicts what
Daniel was saying

Internal functions should *not* check their input conditions.
So you don't want the symmetric handler to check for symmetry
or the 2x2 handler to check sizes.  That's the job of the calling
function and should be tested from there.  It helps ensure the
error reports are correct (they come from the function the user
called, not some internal function they've never heard of)
and removes redundancy in tests.

+1

 

You can test that internal errors arise from the wrong inputs,
but I wouldn't bother.

There are two drawbacks to testing internal functions:

* when the internal functions change without changing the external
  API, then the tests need to change even though it doesn't
  change external public behavior

* requiring tests for internal functions discourages developers
  from breaking natural functional units into tests

I don't get what "breaking natural functional units into tests" means.

 

This would all be mitigated somewhat if the testing framework
were faster.  Right now, with a lot of functions not just
including what they use, a change in a low-level functions requires
a 45 minute unit test run on my end,

Shouldn't we just be testing that the low-level function passes tests? If the low-level function passes tests but higher-level functions fail because of it then the low-level function wasn't tested properly (or it was called for arguments outside its range)

[Avraham Adler]

On Thursday, September 8, 2016 at 11:24:13 AM UTC-4, Krzysztof Sakrejda wrote:
> If the top level function tests that it's branching correctly and the lower-level functions test that they return the correct results then having the top level function test for correct results only duplicates tests. Maybe ideally we'd love to have the top-level functions do that but if the top-level function test fails you have no idea what level the calculation failed at.  We were just trying to hash this out w.r.t. the tests for the Stan refactor and I thought the conclusion was that the function should test it's own code and assume the lower levels are tested separately... 

Faux-programmer here, but I can imagine a case where the lower level function thinks it needs to return something in form A and does so properly, but top level requires B (imaging lower level returning a complex and the upper expecting a real) so it still pays to test that the lower level returns what it is supposed return and the upper level receives what *it* is supposed to receive.

Avi

[Avraham Adler]

On Wednesday, September 7, 2016 at 7:56:23 PM UTC-4, Bob Carpenter wrote:
> We definitely don't want to perform a general eigendecomposition
> to test whether we can apply a more efficient matrix exponential
> algorithm. Is there a fast way to just find the largest eigenvalue?
> If not, you'd want to break the nilpotent case out. I have no
> idea what it's used for and if you don't have an application in
> mind, I'd just recommend skipping it for now.
>
> - Bob
>

Sadly, it appears that the determinant and trace being 0 is necessary but not sufficient condition for nilpotency. The only necessary and sufficient condition I've come across is having all eigenvalues 0.

I'm curious Charles, how are you testing?

[Bob Carpenter]

> On Sep 8, 2016, at 11:24 AM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>
> On Thu, Sep 8, 2016 at 10:55 AM Bob Carpenter <ca...@alias-i.com> wrote:
> -1, since we're voting :-)
>
> The top-level function *must* test [snip] that the branches return the right results.
>
> If the top level function tests that it's branching correctly and the lower-level functions test that they return the correct results then having the top level function test for correct results only duplicates tests.

Exactly my point.

The tests that absolutely must be done are the top-level function
tests for all of its functionality, because those are the functions
that are exposed to the rest of the system. So it's those we need to
make sure keep working when code is optimized or otherwise refactored.

Daniel quite reasonably also wants the lower-level functions tested,
partly so there's some documentation for how they are called, and
partly for a belt-and-suspenders (belt-and-bracers if you're British)
degree of caution.

> Maybe ideally we'd love to have the top-level functions do that but if the top-level function test fails you have no idea what level the calculation failed at.

That's always the case with functions that fail. The functions
don't expose their internals by design, so you have to look
at the code.

The higher-level point is that you want to document a function
thoroughly, including boundary conditions, but not including
internals, and then test vs. the documentation. The problem is
that the encapsulation isn't perfect, so you need to know where
the internal branch points are in the code to make sure you're
testing, say the 2x2 case, because it uses a different algorithm.
That's not going to be part of the interface for the function,
though---the interface just says that it takes a matrix as input
and returns the matrix exponential.

- Bob

[Bob Carpenter]

> On Sep 8, 2016, at 11:28 AM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>
> ...

> * requiring tests for internal functions discourages developers
> from breaking natural functional units into tests
>
>
> I don't get what "breaking natural functional units into tests" means.

Not surprising---it was a typo. I meant breaking the natural
functional units into subroutines.

> This would all be mitigated somewhat if the testing framework
> were faster. Right now, with a lot of functions not just
> including what they use, a change in a low-level functions requires
> a 45 minute unit test run on my end,
>
> Shouldn't we just be testing that the low-level function passes tests? If the low-level function passes tests but higher-level functions fail because of it then the low-level function wasn't tested properly (or it was called for arguments outside its range)

We want to test the higher-level functions thoroughly so that
if their implementation changes, they are still being thoroughly
tested.

- Bob

[Bob Carpenter]

In that particular case, it won't compile because there's
no way to build a real out of a complex. But in general,
the end-to-end test on the top-level function will be testing
that it gets what it needs from the lower-level function.

You can't test the internals of a function because they're
not exposed, so there's no way to get at them. The best you
can do is make sure the top-level tests cover all the internal
logic. That's what I meant by the encapsulation not being
perfect for testing. The tests need to peek into the private
parts of the function to ensure every branch of the function's
logic is being properly tested, including for boundary conditions.

We're nowhere near this level of thorough testing in Stan yet.
We're being better with new functions, but we need to retrofit
a whole bunch of tests to get full test coverage. There are tools
that measure how many of the branches in a function get tested
and we're nowhere near 100%, which is scary.

- Bob

[Bob Carpenter]

> On Sep 8, 2016, at 11:28 AM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>
> ...

> * requiring tests for internal functions discourages developers
> from breaking natural functional units into tests
>
>
> I don't get what "breaking natural functional units into tests" means.

Not surprising---it was a typo. I meant breaking the natural
functional units into subroutines.

> This would all be mitigated somewhat if the testing framework
> were faster. Right now, with a lot of functions not just
> including what they use, a change in a low-level functions requires
> a 45 minute unit test run on my end,
>
> Shouldn't we just be testing that the low-level function passes tests? If the low-level function passes tests but higher-level functions fail because of it then the low-level function wasn't tested properly (or it was called for arguments outside its range)

[Krzysztof Sakrejda]

Ok, thanks for the exposition, it sounds like the top-level tests come first and the low-level tests are nice if we can manage to get them.  Once I finish the test classes for the refactor branch this week I'll try to put some of this down in the Wiki so you don't have to redo the explanation it every time :)  

[Bob Carpenter]

We're going to require the low-level tests because
Daniel wants them.

I'm just trying to explain why we need the high-level
tests and why the low-level functions don't need to
come with bounds checks.

- Bob

[Daniel Lee]

On Thu, Sep 8, 2016 at 11:49 PM, Bob Carpenter <ca...@alias-i.com> wrote:

We're going to require the low-level tests because
Daniel wants them. 

Thanks! I've spent more time than I've wanted in internal functions that have no doc or tests because that's where a bug popped up. It's a lot easier to trap bugs when some tests exists. Even though everyone claims to be around to debug, it's not always true for multiple reasons -- sometimes the symptoms don't seem related but end up in places that aren't tested.

I'm just trying to explain why we need the high-level
tests and why the low-level functions don't need to
come with bounds checks.

Thanks. I'm in complete agreement that low-level functions don't need to check inputs. The boundary conditions should be documented.

Daniel

 

- Bob

> On Sep 8, 2016, at 2:30 PM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>
> Ok, thanks for the exposition, it sounds like the top-level tests come first and the low-level tests are nice if we can manage to get them.  Once I finish the test classes for the refactor branch this week I'll try to put some of this down in the Wiki so you don't have to redo the explanation it every time :)
>
> On Thu, Sep 8, 2016 at 2:05 PM Bob Carpenter <ca...@alias-i.com> wrote:
>
> > On Sep 8, 2016, at 11:28 AM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
> >
> > ...
>
> > * requiring tests for internal functions discourages developers
> > from breaking natural functional units into tests
> >
> >
> > I don't get what "breaking natural functional units into tests" means.
>
> Not surprising---it was a typo.  I meant breaking the natural
> functional units into subroutines.
>
>
> > This would all be mitigated somewhat if the testing framework
> > were faster.  Right now, with a lot of functions not just
> > including what they use, a change in a low-level functions requires
> > a 45 minute unit test run on my end,
> >
> > Shouldn't we just be testing that the low-level function passes tests? If the low-level function passes tests but higher-level functions fail because of it then the low-level function wasn't tested properly (or it was called for arguments outside its range)
>
> We want to test the higher-level functions thoroughly so that
> if their implementation changes, they are still being thoroughly
> tested.
>
> - 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+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.

> 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]

The only necessary and sufficient condition I've come across is having all eigenvalues 0.
I'm curious Charles, how are you testing?

For an NxN matrix A, I computed A *= A, until I either got a zero matrix or I had done the operation (N - 1) times. Not optimal, but easy to code. 

[Avraham Adler]

On Friday, September 9, 2016 at 2:54:45 PM UTC-4, Charles Margossian wrote:
> The only necessary and sufficient condition I've come across is having all eigenvalues 0.
> I'm curious Charles, how are you testing?
>
>
> For an NxN matrix A, I computed A *= A, until I either got a zero matrix or I had done the operation (N - 1) times. Not optimal, but easy to code. 

Actually, I was benchmarking the following code, and the multiple multiplication is not only more than 5 times as fast (for a 10x10 matrix), but more accurate in that it returned TRUE properly for a 10x10 non-obvious nilpotent matrix where the Eigenvalues are subject to floating point issues, I guess, and are not close enough to 0.

Matrix:

> A
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 2 2 2 2 2 2 2 2 2 -9
[2,] 12 1 1 1 1 1 1 1 1 -10
[3,] 1 12 1 1 1 1 1 1 1 -10
[4,] 1 1 12 1 1 1 1 1 1 -10
[5,] 1 1 1 12 1 1 1 1 1 -10
[6,] 1 1 1 1 12 1 1 1 1 -10
[7,] 1 1 1 1 1 12 1 1 1 -10
[8,] 1 1 1 1 1 1 12 1 1 -10
[9,] 1 1 1 1 1 1 1 12 1 -10
[10,] 1 1 1 1 1 1 1 1 12 -10

R code: A <- matrix(c(2, 2, 2, 2, 2, 2, 2, 2, 2, -9, 12, 1, 1, 1, 1, 1, 1, 1, 1, -10, 1, 12, 1, 1, 1, 1, 1, 1, 1, -10, 1, 1, 12, 1, 1, 1, 1, 1, 1, -10, 1, 1, 1, 12, 1, 1, 1, 1, 1, -10, 1, 1, 1, 1, 12, 1, 1, 1, 1, -10, 1, 1, 1, 1, 1, 12, 1, 1, 1, -10, 1, 1, 1, 1, 1, 1, 12, 1, 1, -10, 1, 1, 1, 1, 1, 1, 1, 12, 1, -10, 1, 1, 1, 1, 1, 1, 1, 1, 12, -10), ncol = 10, byrow = TRUE)

C++ Code:

#include <RcppEigen.h>
#include <algorithm>

// [[Rcpp::depends(RcppEigen)]]

using namespace Rcpp;

using Eigen::Map; // 'maps' rather than copies
using Eigen::MatrixXd; // variable size matrix, double precision
using Eigen::VectorXd; // variable size vector, complex double precision
using Eigen::EigenSolver; // general real eigenvalue solvers

const double tol = 1E-12;


// [[Rcpp::export]]
bool Eigentest(Map<MatrixXd> M) {
EigenSolver<MatrixXd> es(M);
VectorXd E = es.eigenvalues().real();
return (E.maxCoeff() < tol);
}

// [[Rcpp::export]]
bool MMULT(MatrixXd M) {
std::size_t n = M.rows();
MatrixXd A = M;
int i = 0;
do {
A *= M;
++i;
} while (A.sum() > tol && i < n);
return(A.sum() < tol);
}

Value results:

> Eigentest(A); MMULT(A)
[1] FALSE
[1] TRUE

Timing Tests:
> microbenchmark(Eigentest(A), MMULT(A), times = 10000, control = list(order = 'block'))
Unit: microseconds
expr min lq mean median uq max neval cld
Eigentest(A) 35.489 36.086 36.563037 36.087 36.385 1377.239 10000 b
MMULT(A) 5.368 5.667 6.356559 5.667 5.965 1324.453 10000 a
>

Thanks,

Avi

[Charles Margossian]

@Avi: Thanks for sharing the code! I'm glad the algorithm is useful -- maybe we'll implement it in Stan after all. 

[Ben Goodrich]

On Friday, September 9, 2016 at 2:54:45 PM UTC-4, Charles Margossian wrote:

The only necessary and sufficient condition I've come across is having all eigenvalues 0.
I'm curious Charles, how are you testing?

For an NxN matrix A, I computed A *= A, until I either got a zero matrix or I had done the operation (N - 1) times. Not optimal, but easy to code. 

Unless we eventually have a parameter type for nilpotent_matrix, I don't see what the point is in worrying about this special case. Either it is a matrix in the parameter block, in which case it will be nilpotent with probability zero or it is a matrix (that should be) in the transformed parameters block, in which case using the Pade approximation once won't hurt us. It is only relevant to optimize the nilpotent case for unoptimized Stan users who repeatedly calculate the expm of a constant nilpotent matrix.

Ben