search

stan-dev / stan

Stan development repository. The master branch contains the current release. The develop branch contains the latest stable development. See the Developer Process Wiki for details.
https://mc-stan.org
BSD 3-Clause "New" or "Revised" License
2.59k stars 369 forks source link

add function mdivide_left_spd and mdivide_right_spd #227

Closed bob-carpenter closed 9 years ago

bob-carpenter commented 11 years ago

mdivide_left_spd is in the code, but not exposed to Stan. Expose it and the right-handed version.

PeterLi2016 commented 10 years ago

Is there a reason these aren't exposed? They have some tests implemented.

bob-carpenter commented 10 years ago

To expose them, we have to decide how to deal with checking that the input truly is symmetric and positive definite. That's probably not necessary for nested uses.

On Jul 30, 2014, at 1:01 PM, Peter Li notifications@github.com wrote:

Is there a reason these aren't exposed? They have some tests implemented.

— Reply to this email directly or view it on GitHub.

syclik commented 9 years ago

I'm suggesting we remove the implementation as part of #1407.

syclik commented 9 years ago

@bgoodri, this is what I need from you:

I need you to give me instances of these matrices:

Are there any other conditions where we should stop A?

Are there any conditions where matrix b should fail?

There are going to be issues just based on inspection. The check_pos_definite comes before checking for a square matrix. I don't know how safe that function is.

bgoodri commented 9 years ago

OK, I can do that. If there are any NaN or +/- Inf, it should fail the check_pos_definite but we may want to have check_pos_definite catch those specifically. This is another case where we want to compute the Eigen::LLT<> and then check it, rather than doing a Cholesky decomposition to see if we can do the same Cholesky decomposition.

On Wed, Mar 25, 2015 at 7:39 PM, Daniel Lee notifications@github.com wrote:

@bgoodri https://github.com/bgoodri, this is what I need from you:

I need you to give me instances of these matrices:

  • matrix A is not symmetric
  • matrix A is not positive definite
  • matrix A is not square
  • matrix A that is ok; matrix b that is not compatible with A Just pick these things out of a hat. I don't care what they are, as long as you provide the actual instances of the matrices.

Are there any other conditions where we should stop A?

  • NaN?
  • Inf? -Inf?
  • zero size?

Are there any conditions where matrix b should fail?

  • NaN?
  • Inf?
  • -Inf?
  • zero size?

There are going to be issues just based on inspection. The check_pos_definite comes before checking for a square matrix. I don't know how safe that function is.

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/stan/issues/227#issuecomment-86253771.

syclik commented 9 years ago

The branch Rob just finished up today has two new function signatures for check_pos_definite: one for LDLT and one for LLT. We can definitely use that.

By the way, it's not safe to assume that our checks catch NaN or that Eigen stops with NaN values.

On Wed, Mar 25, 2015 at 7:43 PM, bgoodri notifications@github.com wrote:

OK, I can do that. If there are any NaN or +/- Inf, it should fail the check_pos_definite but we may want to have check_pos_definite catch those specifically. This is another case where we want to compute the Eigen::LLT<> and then check it, rather than doing a Cholesky decomposition to see if we can do the same Cholesky decomposition.

On Wed, Mar 25, 2015 at 7:39 PM, Daniel Lee notifications@github.com wrote:

@bgoodri https://github.com/bgoodri, this is what I need from you:

I need you to give me instances of these matrices:

  • matrix A is not symmetric
  • matrix A is not positive definite
  • matrix A is not square
  • matrix A that is ok; matrix b that is not compatible with A Just pick these things out of a hat. I don't care what they are, as long as you provide the actual instances of the matrices.

Are there any other conditions where we should stop A?

  • NaN?
  • Inf? -Inf?
  • zero size?

Are there any conditions where matrix b should fail?

  • NaN?
  • Inf?
  • -Inf?
  • zero size?

There are going to be issues just based on inspection. The check_pos_definite comes before checking for a square matrix. I don't know how safe that function is.

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/stan/issues/227#issuecomment-86253771.

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/stan/issues/227#issuecomment-86254307.

bgoodri commented 9 years ago

Eigen doesn't stop but the .diagonal().array() <= 0).any() check should return false if either there are NaNs or divisions by Inf.

On Wed, Mar 25, 2015 at 7:46 PM, Daniel Lee notifications@github.com wrote:

The branch Rob just finished up today has two new function signatures for check_pos_definite: one for LDLT and one for LLT. We can definitely use that.

By the way, it's not safe to assume that our checks catch NaN or that Eigen stops with NaN values.

On Wed, Mar 25, 2015 at 7:43 PM, bgoodri notifications@github.com wrote:

OK, I can do that. If there are any NaN or +/- Inf, it should fail the check_pos_definite but we may want to have check_pos_definite catch those specifically. This is another case where we want to compute the Eigen::LLT<> and then check it, rather than doing a Cholesky decomposition to see if we can do the same Cholesky decomposition.

On Wed, Mar 25, 2015 at 7:39 PM, Daniel Lee notifications@github.com wrote:

@bgoodri https://github.com/bgoodri, this is what I need from you:

I need you to give me instances of these matrices:

  • matrix A is not symmetric
  • matrix A is not positive definite
  • matrix A is not square
  • matrix A that is ok; matrix b that is not compatible with A Just pick these things out of a hat. I don't care what they are, as long as you provide the actual instances of the matrices.

Are there any other conditions where we should stop A?

  • NaN?
  • Inf? -Inf?
  • zero size?

Are there any conditions where matrix b should fail?

  • NaN?
  • Inf?
  • -Inf?
  • zero size?

There are going to be issues just based on inspection. The check_pos_definite comes before checking for a square matrix. I don't know how safe that function is.

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/stan/issues/227#issuecomment-86253771.

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/stan/issues/227#issuecomment-86254307.

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/stan/issues/227#issuecomment-86254856.

syclik commented 9 years ago

@bgoodri, that statement is incorrect. Eigen's transforms swallows both nan and inf. Pop this into src/test/eigen_test.cpp:

#include <gtest/gtest.h>
#include <Eigen/Dense>
#include <limits>

TEST(eigen, ok) {
  Eigen::MatrixXd x(2,2);
  x << 1, 0, 0, 1;

  EXPECT_FALSE((x.diagonal().array() <= 0).any());
}

TEST(eigen, ok2) {
  Eigen::MatrixXd x(2,2);
  x << 0, 0, 0, 0;

  EXPECT_TRUE((x.diagonal().array() <= 0).any());
}

TEST(eigen, nan) {
  Eigen::MatrixXd x(2,2);
  x << 1, 0, 0, std::numeric_limits<double>::quiet_NaN();
  EXPECT_TRUE((x.array() <= 0.0).any());
}

TEST(eigen, nan_ldlt) {
  Eigen::MatrixXd x(2,2);  
  x << 1, 0, 0, std::numeric_limits<double>::quiet_NaN();
  Eigen::LDLT<Eigen::MatrixXd> ldlt = x.ldlt();

  EXPECT_TRUE((ldlt.vectorD().array() <= 0.0).any());
}

TEST(eigen, nan_llt) {
  Eigen::MatrixXd x(2,2);  
  x << 1, 0, 0, std::numeric_limits<double>::quiet_NaN();
  Eigen::LLT<Eigen::MatrixXd> llt = x.llt();

  EXPECT_TRUE((llt.matrixLLT().diagonal().array() <= 0.0).any());
}

TEST(eigen, inf) {
  Eigen::MatrixXd x(2,2);
  x << 1, 0, 0, std::numeric_limits<double>::infinity();
  EXPECT_TRUE((x.array() <= 0.0).any());
}

TEST(eigen, inf_ldlt) {
  Eigen::MatrixXd x(2,2);  
  x << 1, 0, 0, std::numeric_limits<double>::infinity();
  Eigen::LDLT<Eigen::MatrixXd> ldlt = x.ldlt();

  EXPECT_TRUE((ldlt.vectorD().array() <= 0.0).any());
}

TEST(eigen, inf_llt) {
  Eigen::MatrixXd x(2,2);  
  x << 1, 0, 0, std::numeric_limits<double>::infinity();
  Eigen::LLT<Eigen::MatrixXd> llt = x.llt();

  EXPECT_TRUE((llt.matrixLLT().diagonal().array() <= 0.0).any());
}

the output I get:

Running main() from gtest_main.cc
[==========] Running 8 tests from 1 test case.
[----------] Global test environment set-up.
[----------] 8 tests from eigen
[ RUN      ] eigen.ok
[       OK ] eigen.ok (0 ms)
[ RUN      ] eigen.ok2
[       OK ] eigen.ok2 (0 ms)
[ RUN      ] eigen.nan
[       OK ] eigen.nan (0 ms)
[ RUN      ] eigen.nan_ldlt
src/test/eigen_test.cpp:30: Failure
Value of: (ldlt.vectorD().array() <= 0.0).any()
  Actual: false
Expected: true
[  FAILED  ] eigen.nan_ldlt (0 ms)
[ RUN      ] eigen.nan_llt
src/test/eigen_test.cpp:38: Failure
Value of: (llt.matrixLLT().diagonal().array() <= 0.0).any()
  Actual: false
Expected: true
[  FAILED  ] eigen.nan_llt (0 ms)
[ RUN      ] eigen.inf
[       OK ] eigen.inf (0 ms)
[ RUN      ] eigen.inf_ldlt
src/test/eigen_test.cpp:53: Failure
Value of: (ldlt.vectorD().array() <= 0.0).any()
  Actual: false
Expected: true
[  FAILED  ] eigen.inf_ldlt (0 ms)
[ RUN      ] eigen.inf_llt
src/test/eigen_test.cpp:61: Failure
Value of: (llt.matrixLLT().diagonal().array() <= 0.0).any()
  Actual: false
Expected: true
[  FAILED  ] eigen.inf_llt (0 ms)
[----------] 8 tests from eigen (0 ms total)

[----------] Global test environment tear-down
[==========] 8 tests from 1 test case ran. (0 ms total)
[  PASSED  ] 4 tests.
[  FAILED  ] 4 tests, listed below:
[  FAILED  ] eigen.nan_ldlt
[  FAILED  ] eigen.nan_llt
[  FAILED  ] eigen.inf_ldlt
[  FAILED  ] eigen.inf_llt

 4 FAILED TESTS
bgoodri commented 9 years ago

I see what you mean now. We need to change the check from "any non-negative" to "not all positive". Like

!(ldlt.vectorD().array() > 0.0).all();

On Wed, Mar 25, 2015 at 10:33 PM, Daniel Lee notifications@github.com wrote:

@bgoodri https://github.com/bgoodri, that statement is incorrect. Eigen's transforms swallows both nan and inf. Pop this into src/test/eigen_test.cpp:

include <gtest/gtest.h>

include <Eigen/Dense>

include

TEST(eigen, ok) { Eigen::MatrixXd x(2,2); x << 1, 0, 0, 1;

EXPECT_FALSE((x.diagonal().array() <= 0).any()); }

TEST(eigen, ok2) { Eigen::MatrixXd x(2,2); x << 0, 0, 0, 0;

EXPECT_TRUE((x.diagonal().array() <= 0).any()); }

TEST(eigen, nan) { Eigen::MatrixXd x(2,2); x << 1, 0, 0, std::numeric_limits::quiet_NaN(); EXPECT_TRUE((x.array() <= 0.0).any()); }

TEST(eigen, nan_ldlt) { Eigen::MatrixXd x(2,2); x << 1, 0, 0, std::numeric_limits::quiet_NaN(); Eigen::LDLTEigen::MatrixXd ldlt = x.ldlt();

EXPECT_TRUE((ldlt.vectorD().array() <= 0.0).any()); }

TEST(eigen, nan_llt) { Eigen::MatrixXd x(2,2); x << 1, 0, 0, std::numeric_limits::quiet_NaN(); Eigen::LLTEigen::MatrixXd llt = x.llt();

EXPECT_TRUE((llt.matrixLLT().diagonal().array() <= 0.0).any()); }

TEST(eigen, inf) { Eigen::MatrixXd x(2,2); x << 1, 0, 0, std::numeric_limits::infinity(); EXPECT_TRUE((x.array() <= 0.0).any()); }

TEST(eigen, inf_ldlt) { Eigen::MatrixXd x(2,2); x << 1, 0, 0, std::numeric_limits::infinity(); Eigen::LDLTEigen::MatrixXd ldlt = x.ldlt();

EXPECT_TRUE((ldlt.vectorD().array() <= 0.0).any()); }

TEST(eigen, inf_llt) { Eigen::MatrixXd x(2,2); x << 1, 0, 0, std::numeric_limits::infinity(); Eigen::LLTEigen::MatrixXd llt = x.llt();

EXPECT_TRUE((llt.matrixLLT().diagonal().array() <= 0.0).any()); } ``

the output I get:

Running main() from gtest_main.cc [==========] Running 8 tests from 1 test case. [----------] Global test environment set-up. [----------] 8 tests from eigen [ RUN ] eigen.ok [ OK ] eigen.ok (0 ms) [ RUN ] eigen.ok2 [ OK ] eigen.ok2 (0 ms) [ RUN ] eigen.nan [ OK ] eigen.nan (0 ms) [ RUN ] eigen.nan_ldlt src/test/eigen_test.cpp:30: Failure Value of: (ldlt.vectorD().array() <= 0.0).any() Actual: false Expected: true [ FAILED ] eigen.nan_ldlt (0 ms) [ RUN ] eigen.nan_llt src/test/eigen_test.cpp:38: Failure Value of: (llt.matrixLLT().diagonal().array() <= 0.0).any() Actual: false Expected: true [ FAILED ] eigen.nan_llt (0 ms) [ RUN ] eigen.inf [ OK ] eigen.inf (0 ms) [ RUN ] eigen.inf_ldlt src/test/eigen_test.cpp:53: Failure Value of: (ldlt.vectorD().array() <= 0.0).any() Actual: false Expected: true [ FAILED ] eigen.inf_ldlt (0 ms) [ RUN ] eigen.inf_llt src/test/eigen_test.cpp:61: Failure Value of: (llt.matrixLLT().diagonal().array() <= 0.0).any() Actual: false Expected: true [ FAILED ] eigen.inf_llt (0 ms) [----------] 8 tests from eigen (0 ms total)

[----------] Global test environment tear-down [==========] 8 tests from 1 test case ran. (0 ms total) [ PASSED ] 4 tests. [ FAILED ] 4 tests, listed below: [ FAILED ] eigen.nan_ldlt [ FAILED ] eigen.nan_llt [ FAILED ] eigen.inf_ldlt [ FAILED ] eigen.inf_llt

4 FAILED TESTS

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/stan/issues/227#issuecomment-86302944.

syclik commented 9 years ago

This issue was moved to stan-dev/math#22