adding input checks for cholesky factors #3002

syclik · 2024-01-13T01:15:16Z

Summary

We weren't checking cholesky factors in the code in multivariate distributions. This appeared in a user-facing RNG generating NaN for all quantities.

The fix here is to validate input arguments.

Tests

No new tests.

For anyone reviewing, if it makes sense, please let me know what should be tested.

Side Effects

Some existing models may no longer work if they accepted matrixes and computed something (even if it's non-sensical).

Release notes

Added checks for Cholesky factors in multi-variate distributions.

Checklist

Copyright holder: Daniel Lee

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

spinkney · 2024-01-17T01:19:27Z

@syclik I fixed the issue and pushed. Also I updated the tests to reflect the cholesky check and fixed a ref to nu.

spinkney

The fix was to switch the factors and use mdivide_right_tri_low. I misread the paper as having the math on a lower tri factor.

I also updated the test to check that it gives, on average, the same values as inv_wishart_rng for the full covariance matrix.

syclik · 2024-01-17T01:24:54Z

Thank you!!! I was going to get around to it soon. Seth replied on slack. Sounded like there was another issue that’s different from input checking in Julia.

…

On Tue, Jan 16, 2024 at 8:19 PM Sean Pinkney ***@***.***> wrote: @syclik <https://github.com/syclik> I fixed the issue and pushed. Also I updated the tests to reflect the cholesky check and fixed a ref to nu. — Reply to this email directly, view it on GitHub <#3002 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AADH6F2ICJA57QDX652DXV3YO4RKVAVCNFSM6AAAAABBY5PQK2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQOJUG43TSOJQHA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

syclik · 2024-01-17T04:46:27Z

test/unit/math/opencl/rev/multi_normal_cholesky_lpdf_test.cpp

@@ -105,7 +105,7 @@ TEST(ProbDistributionsMultiNormalCholesky, opencl_matches_cpu_small) {
  std::vector<Eigen::VectorXd> mu3{mu1, mu2};
  std::vector<Eigen::RowVectorXd> mu4{mu2, mu1};
  Eigen::MatrixXd L(N, N);
-  L << 1, 0, 0, 2, 3, 0, -4, 5, -6;
+  L << 4, 0, 0, -1, 1, 0, -4, -2, 3;


Note: the original code is not a valid cholesky factor! The test was passing incorrectly.

spinkney · 2024-01-17T10:36:35Z

@syclik the test test/unit/math/mix/prob_1_test is failing for inv_wishart_cholesky_lpdf on ProbDistributionsInvWishartCholesky.matvar. Is there a way for us to loosen the tolerances or fix this test?

unknown file: Failure

C++ exception with description "inv_wishart_cholesky_lpdf: Cholesky random variable is not lower triangular; Cholesky random variable[1,2]=1.81664e-05" thrown in the test body.

syclik · 2024-01-17T11:28:41Z

That's a pretty big value. Let me see what's going on with that test.

spinkney · 2024-01-17T11:31:44Z

The issue is with this test

  Eigen::MatrixXd y22(2, 2);
  y22 << 1.0, 0.1, 0.1, 1.5;
  Eigen::MatrixXd L_Y22 = y22.llt().matrixL();
  Eigen::MatrixXd Sigma22(2, 2);
  Sigma22 << 1.5, 0.5, 0.5, 1.1;
  Eigen::MatrixXd L_S22 = Sigma22.llt().matrixL();
  stan::test::expect_ad(f, L_Y22, dof, L_S22);

Finite differences incorrectly perturbs the upper triangular element of the 2 x 2 matrix. Printing out the values of both the cholesky factor matrices shows that the inputs are lower triangular and the 1, 2 element is 0.

This error will again happen with the wishart_cholesky_lpdf in prob_3_test.cpp.

@WardBrian @SteveBronder how do we update these tests for cholesky factors? Or should we just remove the check (in both the inv and regular wishart cholesky lpdfs) which is what @syclik is trying to add here?

andrjohns · 2024-01-17T12:04:00Z

You could change the test functor to call cholesky_factor_constrain on the input before it calls the lpdf

syclik · 2024-01-17T12:05:18Z

I was looking for an equivalent expect_ad_ that accounted for cholesky factors. I can imagine either building a new one, updating f so that it passes if the term is not a cholesky factor, or passing in the matrix and using f to create cholesky factors of the matrix (I think this may still have weird errors if the matrix does not have good properties). I’ll search for how we deal with other constrained types (like symmetric matrix).

…

On Wed, Jan 17, 2024 at 6:31 AM Sean Pinkney ***@***.***> wrote: The issue is with this test Eigen::MatrixXd y22(2, 2); y22 << 1.0, 0.1, 0.1, 1.5; Eigen::MatrixXd L_Y22 = y22.llt().matrixL(); Eigen::MatrixXd Sigma22(2, 2); Sigma22 << 1.5, 0.5, 0.5, 1.1; Eigen::MatrixXd L_S22 = Sigma22.llt().matrixL(); stan::test::expect_ad(f, L_Y22, dof, L_S22); Finite differences incorrectly perturbs the upper triangular element of the 2 x 2 matrix. Printing out the values of both the cholesky factor matrices shows that the inputs are lower triangular and the 1, 2 element is 0. This error will again happen with the wishart_cholesky_lpdf in prob_3_test.cpp. @WardBrian <https://github.com/WardBrian> @SteveBronder <https://github.com/SteveBronder> how do we update these tests for cholesky factors? Or should we just remove the check (in both the inv and regular wishart cholesky lpdfs) which is what @syclik <https://github.com/syclik> is trying to add here? — Reply to this email directly, view it on GitHub <#3002 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AADH6F2V2MBTFMM2CEX47Q3YO6ZCXAVCNFSM6AAAAABBY5PQK2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQOJVGYZDENRSGQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

syclik · 2024-01-17T14:20:23Z

Update: ignore this message. I was incorrect.

@spinkney, just FYI, it's not where you thought. I added prints (because we can't tell where it failed from the output). It's in the first call to expect_ad_matvar()

stan::test::expect_ad_matvar(f, L_Y22, dof, L_S22);

syclik · 2024-01-17T14:27:54Z

@spinkney, ignore the last message; it fails before that. I misread my prints and thought it was further along.

syclik · 2024-01-17T14:35:50Z

I tried writing a new functor and the test fails to compile.

  auto f2 = [](const auto& L_Y, const auto& dof, const auto& Sigma) {
    auto L_S = Sigma.llt().matrixL();
    return stan::math::inv_wishart_cholesky_lpdf(L_Y, dof, L_S);
  };
...
stan::test::expect_ad(f2, L_Y22, dof, Sigma22);

Results in a compiler novel with the first error message:

In file included from test/unit/math/mix/prob/inv_wishart_cholesky_test.cpp:1:
In file included from ./test/unit/math/test_ad.hpp:4:
In file included from ./stan/math/mix.hpp:4:
In file included from ./stan/math/mix/meta.hpp:6:
In file included from ./stan/math/fwd/core.hpp:4:
In file included from ./stan/math/fwd/core/fvar.hpp:4:
In file included from ./stan/math/prim/meta.hpp:72:
In file included from ./stan/math/prim/meta/append_return_type.hpp:4:
In file included from ./stan/math/prim/fun/Eigen.hpp:22:
In file included from lib/eigen_3.4.0/Eigen/Dense:1:
In file included from lib/eigen_3.4.0/Eigen/Core:309:
lib/eigen_3.4.0/Eigen/src/Core/Ref.h:35:82: error: no member named 'InnerStrideAtCompileTime' in 'Eigen::TriangularView<const Eigen::Matrix<double, -1, -1>, 1>'
                      || int(StrideType::InnerStrideAtCompileTime)==int(Derived::InnerStrideAtCompileTime)
                                                                        ~~~~~~~~~^
./stan/math/prim/meta/ref_type.hpp:43:15: note: in instantiation of template class 'Eigen::internal::traits<Eigen::Ref<Eigen::Matrix<double, -1, -1>>>::match<Eigen::TriangularView<const Eigen::Matrix<double, -1, -1>, 1>>' requested here
              template match<T_dec>::MatchAtCompileTime
              ^
./stan/math/prim/meta/ref_type.hpp:55:1: note: in instantiation of template class 'stan::ref_type_if<true, Eigen::TriangularView<const Eigen::Matrix<double, -1, -1>, 1>>' requested here
using ref_type_t = typename ref_type_if<true, T>::type;
^
./stan/math/prim/prob/inv_wishart_cholesky_lpdf.hpp:48:21: note: in instantiation of template type alias 'ref_type_t' requested here
  using T_L_S_ref = ref_type_t<T_scale>;
                    ^
./stan/math/prim/prob/inv_wishart_cholesky_lpdf.hpp:91:10: note: in instantiation of function template specialization 'stan::math::inv_wishart_cholesky_lpdf<false, Eigen::Matrix<double, -1, -1>, double, Eigen::TriangularView<const Eigen::Matrix<double, -1, -1>, 1>, nullptr, nullptr>' requested here
  return inv_wishart_cholesky_lpdf<false>(L_Y, nu, L_S);
         ^
test/unit/math/mix/prob/inv_wishart_cholesky_test.cpp:9:24: note: in instantiation of function template specialization 'stan::math::inv_wishart_cholesky_lpdf<Eigen::Matrix<double, -1, -1>, double, Eigen::TriangularView<const Eigen::Matrix<double, -1, -1>, 1>>' requested here
    return stan::math::inv_wishart_cholesky_lpdf(L_Y, dof, L_S);
                       ^
./test/unit/math/test_ad.hpp:701:13: note: in instantiation of function template specialization 'stan::test::internal::expect_ad_helper<(lambda at test/unit/math/mix/prob/inv_wishart_cholesky_test.cpp:7:13), (lambda at ./test/unit/math/test_ad.hpp:696:13), Eigen::Matrix<double, -1, -1>, double, Eigen::Matrix<double, -1, -1>>' requested here
  internal::expect_ad_helper(tols, f, g1, serialize_args(x1), x1, x2, x3);
            ^
./test/unit/math/test_ad.hpp:1232:13: note: in instantiation of function template specialization 'stan::test::internal::expect_ad_vvv<(lambda at test/unit/math/mix/prob/inv_wishart_cholesky_test.cpp:7:13), Eigen::Matrix<double, -1, -1>, double, Eigen::Matrix<double, -1, -1>>' requested here
  internal::expect_ad_vvv(tols, f, x1, x2, x3);
            ^
./test/unit/math/test_ad.hpp:1252:3: note: in instantiation of function template specialization 'stan::test::expect_ad<(lambda at test/unit/math/mix/prob/inv_wishart_cholesky_test.cpp:7:13), Eigen::Matrix<double, -1, -1>, double, Eigen::Matrix<double, -1, -1>>' requested here
  expect_ad(tols, f, x1, x2, x3);
  ^
test/unit/math/mix/prob/inv_wishart_cholesky_test.cpp:36:15: note: in instantiation of function template specialization 'stan::test::expect_ad<(lambda at test/unit/math/mix/prob/inv_wishart_cholesky_test.cpp:7:13), Eigen::Matrix<double, -1, -1>, double, Eigen::Matrix<double, -1, -1>>' requested here
  stan::test::expect_ad(f2, L_Y22, dof, Sigma22);
              ^

WardBrian · 2024-01-17T14:57:44Z

@spinkney I think what @andrjohns suggested is the best idea, treat the function f as a composition of a constraint and the function you want to test.

This is what we do for e.g. eigenvalues_sym:

  auto f = [](const auto& y) {
    // need to maintain symmetry for finite diffs
    auto a = ((y + y.transpose()) * 0.5).eval();
    return stan::math::eigenvalues_sym(a);
  };

spinkney · 2024-01-17T15:50:33Z

@WardBrian I don't know how to update the generated jumbo test to have this new functor automatically included in the tests.

Secondly, using Eigen's llt() results in a ton of template errors (see @syclik comment). If I use stan::math::cholesky_decompose it gives a non-square matrix (I don't understand that).

WardBrian · 2024-01-17T15:54:17Z

The jumbo tests are literally just compiling multiple existing tests at once to save on time. If you update the individual unit test files, the jumbos will be updated.

The only place we're actually generating test code is in the expression tests. If those also need customization, you can do that by editing the dictionaries in https://github.com/stan-dev/math/blob/develop/test/sig_utils.py

WardBrian · 2024-01-17T16:09:50Z

@syclik's code can be made to compile if you force evaluation of the expression, like so:

  auto f = [](const auto& L_Y, const auto& dof, const auto& Sigma) {
    auto L_S = Sigma.llt().matrixL();
    stan::plain_type_t<decltype(L_S)> L_S_plain = L_S;
    return stan::math::inv_wishart_cholesky_lpdf(L_Y, dof, L_S_plain);
  };

To me this seems like something we'd want the lpdf to be able to handle itself...

However, you would also need to comment out the expect_ad_matvar tests, because we don't define .llt() on var<matrix>s. We do have a specialization of cholesky_decompose for them

WardBrian · 2024-01-17T16:14:19Z

Using this as f compiles and passes for me:

  auto f = [](const auto& L_Y, const auto& dof, const auto& Sigma) {
    auto L_S = stan::math::cholesky_decompose(Sigma);
    stan::plain_type_t<decltype(L_S)> L_S_plain = L_S;
    return stan::math::inv_wishart_cholesky_lpdf(L_Y, dof, L_S);
  };

We'd probably also want to use something like stan::test::square_test_matrices to expand this testing now that it would work on general square matrices

spinkney · 2024-01-17T21:51:22Z

Using this as f compiles and passes for me:

  auto f = [](const auto& L_Y, const auto& dof, const auto& Sigma) {
    auto L_S = stan::math::cholesky_decompose(Sigma);
    stan::plain_type_t<decltype(L_S)> L_S_plain = L_S;
    return stan::math::inv_wishart_cholesky_lpdf(L_Y, dof, L_S);
  };

Is L_S in the last line supposed to be L_S_plain? I tried both and I still experience test failures.

syclik · 2024-01-19T16:17:18Z

I'm trying to work through the expression tests that failed.

WardBrian · 2024-01-19T19:49:20Z

To enable this to be in 4.8.1 without a ton of extra hassle I've created a branch called staging/v4.8.1 and changed the target of the PR to that one. Unfortunately this branch had already had some things from develop merged in, so if you want this as a candidate for inclusion those would need to be removed from this branch's history.

If not, you can just change the target branch back to develop and we will delete the staging/v4.8.1 branch

…g functions

…. New code properly catches the error.

syclik · 2024-01-19T21:08:04Z

Closing. We can use the other PR #3007.

syclik · 2024-01-19T21:17:37Z

Please do not merge this PR. I reopened to kick off tests. #3007 isn't running tests.

syclik · 2024-01-20T18:13:44Z

Closing this PR. #3007 is the PR that should be merged.

syclik requested a review from spinkney January 13, 2024 01:15

spinkney previously approved these changes Jan 17, 2024

View reviewed changes

spinkney dismissed their stale review via a3b3e6b January 17, 2024 01:27

spinkney previously approved these changes Jan 17, 2024

View reviewed changes

spinkney mentioned this pull request Jan 17, 2024

Release/v4.8.0 #3004

Closed

4 tasks

syclik dismissed spinkney’s stale review via 6b767a0 January 17, 2024 04:05

syclik commented Jan 17, 2024

View reviewed changes

serban-nicusor-toptal mentioned this pull request Jan 17, 2024

Release CmdStan 2.34 stan-dev/cmdstan#1230

Closed

26 tasks

syclik force-pushed the feature/check_cholesky_factor branch 3 times, most recently from eba3fc3 to 680b1a9 Compare January 18, 2024 11:55

WardBrian changed the base branch from develop to staging/v4.8.1 January 19, 2024 19:47

syclik and others added 4 commits January 19, 2024 16:00

adding checks for cholesky factors; we found problems with user-facin…

66fff6d

…g functions

Update wishart_cholesky_rng_test.cpp

994e6bd

updating test; test was incorrect. It passed in a non-cholesky factor…

eb5d9b5

…. New code properly catches the error.

passing expression tests

ea31dad

syclik force-pushed the feature/check_cholesky_factor branch from 2565ff0 to ea31dad Compare January 19, 2024 21:01

syclik mentioned this pull request Jan 19, 2024

adding input checks for cholesky factors #3007

Merged

4 tasks

syclik closed this Jan 19, 2024

syclik reopened this Jan 19, 2024

syclik changed the base branch from staging/v4.8.1 to develop January 19, 2024 21:16

syclik and others added 2 commits January 19, 2024 22:20

passing expression tests

40022fc

[Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1

f8611c8

syclik closed this Jan 20, 2024

syclik deleted the feature/check_cholesky_factor branch January 20, 2024 18:14

WardBrian mentioned this pull request Jan 24, 2024

Release/v4.8.1 #3011

Merged

4 tasks

syclik mentioned this pull request Feb 6, 2024

updating opencl multi_normal_cholesky_lpdf to match the other implementation #3021

Open

4 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

adding input checks for cholesky factors #3002

adding input checks for cholesky factors #3002

syclik commented Jan 13, 2024 •

edited by spinkney

Loading

spinkney commented Jan 17, 2024

spinkney left a comment

syclik commented Jan 17, 2024 via email

syclik Jan 17, 2024

spinkney commented Jan 17, 2024

syclik commented Jan 17, 2024

spinkney commented Jan 17, 2024

andrjohns commented Jan 17, 2024

syclik commented Jan 17, 2024 via email

syclik commented Jan 17, 2024 •

edited

Loading

syclik commented Jan 17, 2024

syclik commented Jan 17, 2024

WardBrian commented Jan 17, 2024

spinkney commented Jan 17, 2024

WardBrian commented Jan 17, 2024

WardBrian commented Jan 17, 2024

WardBrian commented Jan 17, 2024

spinkney commented Jan 17, 2024

syclik commented Jan 19, 2024

WardBrian commented Jan 19, 2024

syclik commented Jan 19, 2024

syclik commented Jan 19, 2024

syclik commented Jan 20, 2024

adding input checks for cholesky factors #3002

adding input checks for cholesky factors #3002

Conversation

syclik commented Jan 13, 2024 • edited by spinkney Loading

Summary

Tests

Side Effects

Release notes

Checklist

spinkney commented Jan 17, 2024

spinkney left a comment

Choose a reason for hiding this comment

syclik commented Jan 17, 2024 via email

syclik Jan 17, 2024

Choose a reason for hiding this comment

spinkney commented Jan 17, 2024

syclik commented Jan 17, 2024

spinkney commented Jan 17, 2024

andrjohns commented Jan 17, 2024

syclik commented Jan 17, 2024 via email

syclik commented Jan 17, 2024 • edited Loading

syclik commented Jan 17, 2024

syclik commented Jan 17, 2024

WardBrian commented Jan 17, 2024

spinkney commented Jan 17, 2024

WardBrian commented Jan 17, 2024

WardBrian commented Jan 17, 2024

WardBrian commented Jan 17, 2024

spinkney commented Jan 17, 2024

syclik commented Jan 19, 2024

WardBrian commented Jan 19, 2024

syclik commented Jan 19, 2024

syclik commented Jan 19, 2024

syclik commented Jan 20, 2024

syclik commented Jan 13, 2024 •

edited by spinkney

Loading

syclik commented Jan 17, 2024 •

edited

Loading