Feature/matrix input to mvn [WIP]

[adamhaber]

Summary

Allow multivariate normal distribution to take in matrices for y and mu arguments (see #2532)

Tests

Added a simple initial test for lpdf and rng, will add others (matching the tests for the non-matrix mvn?) if the current implementation makes sense.

Side Effects

None.

Release notes

Allow multivariate normal distribution to take in matrices for y and mu arguments.

Checklist

• Math issue Allow multivariate distributions to take in matrices #2532

• Copyright holder: Adam Haber

  The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
  - Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
  - Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

• the basic tests are passing

  • unit tests pass (to run, use: ./runTests.py test/unit)
  • header checks pass, (make test-headers)
  • dependencies checks pass, (make test-math-dependencies)
  • docs build, (make doxygen)
  • code passes the built in C++ standards checks (make cpplint)

• the code is written in idiomatic C++ and changes are documented in the doxygen

• the new changes are tested

[andrjohns]

For the _lpdf functions, would it be easier to just call to_vector on the inputs? Then the existing vector inputs are unchanged but the matrix inputs are mapped to vectors and the rest of the function proceeds as normal? That might allow for more to be vectorised?

[adamhaber]

Do you mean transforming y and mu from [N,M] matrices to N*M vectors?

[andrjohns]

I was initially, but now I realise that doesn't make much sense in this context, where the y/mu dimensions need to match sigma's. What i was thinking would only really make sense for passing matrices to the univariate distributions

[stan-buildbot]

---

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
gp_pois_regr/gp_pois_regr.stan	3.1	3.07	1.01	0.85% faster
low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.99	-0.96% slower
eight_schools/eight_schools.stan	0.11	0.12	0.98	-1.66% slower
gp_regr/gp_regr.stan	0.16	0.16	0.99	-0.86% slower
irt_2pl/irt_2pl.stan	5.92	5.91	1.0	0.25% faster
performance.compilation	87.93	86.62	1.02	1.5% faster
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	8.57	8.61	0.99	-0.52% slower
pkpd/one_comp_mm_elim_abs.stan	29.31	29.49	0.99	-0.61% slower
sir/sir.stan	131.72	132.27	1.0	-0.42% slower
gp_regr/gen_gp_data.stan	0.03	0.03	0.99	-1.21% slower
low_dim_gauss_mix/low_dim_gauss_mix.stan	3.0	3.0	1.0	0.06% faster
pkpd/sim_one_comp_mm_elim_abs.stan	0.39	0.39	1.02	2.19% faster
arK/arK.stan	2.56	2.53	1.01	1.17% faster
arma/arma.stan	0.65	0.65	1.01	1.32% faster
garch/garch.stan	0.64	0.64	1.0	0.46% faster
Mean result: 1.00115771047

---

Jenkins Console Log
Blue Ocean
Commit hash: 03de55a

---

Machine information

ProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010

CPU:
Intel(R) Xeon(R) CPU E5-1680 v2 @ 3.00GHz

G++:
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

Clang:
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

[SteveBronder]

Thanks Adam! I made an example of what I think would be nice here you can check out in the diff below

What we want to do is utilize the matrix inputs for more efficient code than the standard mvn. So in the below instead of having to call trace_inv_quad_form_ldlt for each row we can just call it once. Would you be able to take something like I have in the below and apply it to the other multivariate normal distributions? Ideally this would be rad to have for all the multivariate distributions but we can do these piece by piece

https://github.com/stan-dev/math/compare/review/mvn-matrix?expand=1

[syclik]

Closing this for now since it's still in draft. @adamhaber, please reopen if it's active.

[syclik]

Sorry... I didn't realize it wasn't a draft. Reopening.

[syclik]

Thanks for the PR!

Please add a unit test for different sizes. I believe it should throw an std::illegal_argument exception so that sampling stops.

[syclik]

Could you take away the using lp_type... statement on 111 and declare lp as that type? It's easier to read when that's together.

[syclik]

Is this code safe? What happens if y and mu have different lengths? Won't this just read off the end?

I'd like to see this fixed before it gets into the code base. Perhaps use a check_matching_sizes or something similar?

[syclik]

Please add a test where y and mu have different number of rows. This will trigger a failure and we should know what that is through a unit test.

[syclik]

Sorry, it is a draft! I'm closing it until it's ready to be reopened.