add function mdivide_left_spd and mdivide_right_spd

[syclik]

From @bob-carpenter on September 20, 2013 20:14

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

• mdivide_right_spd, defined via transpose
• manual doc
• expose function_signatures
• unit tests
• test models

Copied from original issue: stan-dev/stan#227

[syclik]

From @PeterLi2016 on July 30, 2014 17:1

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

[syclik]

From @bob-carpenter on July 30, 2014 17:14

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.

• Bob

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]

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

[syclik]

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

[syclik]

From @bgoodri on March 25, 2015 23:43

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
stan-dev/stan#227 (comment).

[syclik]

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
stan-dev/stan#227 (comment).

—
Reply to this email directly or view it on GitHub
stan-dev/stan#227 (comment).

[syclik]

From @bgoodri on March 25, 2015 23:53

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
stan-dev/stan#227 (comment).

—
Reply to this email directly or view it on GitHub
stan-dev/stan#227 (comment).

—
Reply to this email directly or view it on GitHub
stan-dev/stan#227 (comment).

[syclik]

@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

[syclik]

From @bgoodri on March 26, 2015 3:25

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
stan-dev/stan#227 (comment).

[syclik]

@bgoodri, want to add these tests so we can get these functions exposed?