Skip to content

Commit

Permalink
clean up matrix inverse for fvar
Browse files Browse the repository at this point in the history
  • Loading branch information
@PeterLi2016
PeterLi2016 committed May 24, 2013
1 parent 0c87b34 commit 2b0107e
Show file tree
Hide file tree
Showing 3 changed files with 99 additions and 14 deletions.
7 changes: 6 additions & 1 deletion src/stan/agrad/fwd/matrix/inverse.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,8 @@
#include <stan/agrad/fwd/operator_multiplication.hpp>
#include <stan/math/matrix/multiply.hpp>
#include <stan/agrad/fwd/matrix/to_fvar.hpp>
#include <stan/agrad/fwd/matrix/multiply.hpp>
#include <stan/math/matrix/inverse.hpp>

namespace stan {
namespace agrad {
Expand All @@ -18,6 +20,8 @@ namespace stan {
Eigen::Matrix<fvar<T>,R,C>
inverse(const Eigen::Matrix<fvar<T>, R, C>& m) {
using stan::math::multiply;
using stan::agrad::multiply;
using stan::math::inverse;
stan::math::validate_square(m, "inverse");
Eigen::Matrix<T,R,C> m_deriv(m.rows(), m.cols());
Eigen::Matrix<T,R,C> m_inv(m.rows(), m.cols());
Expand All @@ -29,7 +33,8 @@ namespace stan {
}
}

m_inv = m_inv.inverse();
m_inv = inverse(m_inv);

m_deriv = -1 * multiply(multiply(m_inv, m_deriv), m_inv);

return to_fvar(m_inv, m_deriv);
Expand Down
99 changes: 93 additions & 6 deletions src/stan/agrad/fwd/matrix/mdivide_left.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,10 @@
#include <stan/agrad/fwd/matrix/typedefs.hpp>
#include <stan/agrad/fwd/matrix/inverse.hpp>
#include <stan/agrad/fwd/matrix/multiply.hpp>
#include <stan/math/matrix/multiply.hpp>
#include <stan/agrad/fwd/matrix/to_fvar.hpp>
#include <stan/math/matrix/inverse.hpp>
#include <stan/agrad/fwd/matrix/inverse.hpp>

namespace stan {
namespace agrad {
Expand All @@ -20,10 +24,43 @@ namespace stan {
mdivide_left(const Eigen::Matrix<fvar<T1>,R1,C1> &A,
const Eigen::Matrix<fvar<T2>,R2,C2> &b) {

using stan::agrad::multiply;
using stan::math::multiply;
using stan::agrad::inverse;
using stan::math::inverse;
stan::math::validate_square(A,"mdivide_left");
stan::math::validate_multiplicable(A,b,"mdivide_left");

return stan::agrad::multiply(stan::agrad::inverse(A), b);

Eigen::Matrix<fvar<T1>,R1,C1> fvar_inv_A;
fvar_inv_A = stan::agrad::inverse(A);

This comment has been minimized.

Sign in to view

Copy link
@bob-carpenter

bob-carpenter May 29, 2013

Contributor

We should not be taking inverses! A major application of left division like this is to avoid taking matrix inverses.

This code needs to be rewritten to follow the way the reverse-mode code works in terms of the linear algebra solvers.


Eigen::Matrix<typename stan::return_type<T1,T2>,R1,C2> inv_A_mult_b;
Eigen::Matrix<typename stan::return_type<T1,T2>,R1,C2> val;
Eigen::Matrix<T1,R1,C1> val_A;
Eigen::Matrix<T1,R1,C2> inv_A;
Eigen::Matrix<T1,R1,C1> deriv_A;
Eigen::Matrix<T2,R2,C2> val_b;
Eigen::Matrix<T2,R2,C2> deriv_b;

for (int i = 0; i < A.rows(); i++) {
for(int j = 0; j < A.cols(); j++) {
inv_A(i,j) = fvar_inv_A(i,j).val_;
val_A(i,j) = A(i,j).d_;
deriv_A(i,j) = A(i,j).d_;
val_b(i,j) = b(i,j).val_;
deriv_b(i,j) = b(i,j).d_;
}
}

inv_A_mult_b = multiply(inverse(val_A), val_b);

Eigen::Matrix<typename stan::return_type<T1,T2>,R1,C2> deriv;
deriv = multiply(inv_A, val_b);
deriv = multiply(deriv_A, deriv);
deriv = multiply(inv_A, deriv_b - deriv);

return stan::agrad::to_fvar(inv_A_mult_b, deriv);
}

template <typename T1, typename T2, int R1,int C1,int R2,int C2>
Expand All @@ -32,10 +69,38 @@ namespace stan {
mdivide_left(const Eigen::Matrix<fvar<T1>,R1,C1> &A,
const Eigen::Matrix<T2,R2,C2> &b) {

using stan::agrad::multiply;
using stan::math::multiply;
using stan::agrad::inverse;
using stan::math::inverse;
stan::math::validate_square(A,"mdivide_left");
stan::math::validate_multiplicable(A,b,"mdivide_left");

return stan::agrad::multiply(stan::agrad::inverse(A), b);

Eigen::Matrix<fvar<T1>,R1,C1> fvar_inv_A;
fvar_inv_A = stan::agrad::inverse(A);

Eigen::Matrix<typename stan::return_type<T1,T2>,R1,C2> inv_A_mult_b;
Eigen::Matrix<typename stan::return_type<T1,T2>,R1,C2> val;
Eigen::Matrix<T1,R1,C1> val_A;
Eigen::Matrix<T1,R1,C1> deriv_A;
Eigen::Matrix<T1,R1,C1> inv_A;

for (int i = 0; i < A.rows(); i++) {
for(int j = 0; j < A.cols(); j++) {

This comment has been minimized.

Sign in to view

Copy link
@randommm

randommm May 29, 2013

Member

I think the correct order is the big loop i < A.cols(); and the smaller j < A.rows();.

Anyway, I found using xs.data() is noticeable quicker for big matrices http://stackoverflow.com/questions/16283000/most-efficient-way-to-loop-thought-eigen-matrix/16287999

So, I think here, we could have something like:

fvar<T1> * ABegin = A.data();
fvar<T1> * fvar_inv_ABegin = fvar_inv_A.data();
T1 * inv_ABegin = fvar_inv_A.data();
T1 * val_ABegin = fvar_inv_A.data();
T1 * deriv_ABegin = fvar_inv_A.data();
for (size_t i = 0, size = A.size(); i < size; i++)
{
   inv_ABegin[i] = fvar_inv_ABegin[i].val_;
   val_ABegin[i] = ABegin[i].d_;
   deriv_ABegin[i] = ABegin[i].d_;
}

Also, at least for me, compiling with -Wall and use int instead of size_t in the code there, you'll get a warning: comparison between signed and unsigned integer expressions [-Wsign-compare]. Jenkins received such warning for sort functions I made using size_t, so maybe the behavior is different on Windows, so I don't know about this...

This comment has been minimized.

Sign in to view

Copy link
@bob-carpenter

bob-carpenter May 29, 2013

Contributor
  1. Yes, we definitely want to loop in column-major order, because we're using the default Eigen ordering.
  2. We've been pretty sloppy on indexing, using int for Eigen and size_t for vector. But we're trying to do better. For Eigen indexing, you should use
Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic>::Index

instead of size_t or int.

For standard vector indexing, you should use

std::vector<T>::size_type

where T is the type of the vector.
inv_A(i,j) = fvar_inv_A(i,j).val_;
val_A(i,j) = A(i,j).d_;
deriv_A(i,j) = A(i,j).d_;
}
}

inv_A_mult_b = multiply(inverse(val_A), b);

Eigen::Matrix<typename stan::return_type<T1,T2>,R1,C2> deriv;
deriv = multiply(inv_A, b);
deriv = multiply(deriv_A, deriv);
deriv = -multiply(inv_A, deriv);

return stan::agrad::to_fvar(inv_A_mult_b, deriv);
}

template <typename T1, typename T2, int R1,int C1,int R2,int C2>
Expand All @@ -44,10 +109,32 @@ namespace stan {
mdivide_left(const Eigen::Matrix<T1,R1,C1> &A,
const Eigen::Matrix<fvar<T2>,R2,C2> &b) {

using stan::agrad::multiply;
using stan::math::multiply;
using stan::agrad::inverse;
using stan::math::inverse;
stan::math::validate_square(A,"mdivide_left");
stan::math::validate_multiplicable(A,b,"mdivide_left");

return stan::agrad::multiply(stan::agrad::inverse(to_fvar(A)), b);

Eigen::Matrix<T1,R1,C1> inv_A;
inv_A = inverse(A);

Eigen::Matrix<typename stan::return_type<T1,T2>,R1,C2> inv_A_mult_b;
Eigen::Matrix<T2,R2,C2> deriv_b;
Eigen::Matrix<T2,R2,C2> val_b;

for (int i = 0; i < A.rows(); i++) {
for(int j = 0; j < A.cols(); j++) {
deriv_b(i,j) = b(i,j).d_;
val_b(i,j) = b(i,j).val_;
}
}

inv_A_mult_b = multiply(inv_A, val_b);

Eigen::Matrix<typename stan::return_type<T1,T2>,R1,C2> deriv;
deriv = multiply(inv_A, deriv_b);
return stan::agrad::to_fvar(inv_A_mult_b, deriv);
}
}
}
Expand Down
7 changes: 0 additions & 7 deletions src/stan/agrad/fwd/matrix/multiply.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,13 +14,6 @@

namespace stan {
namespace agrad {

template <typename T1, typename T2>
inline
typename stan::return_type<T1,T2>::type
multiply(const T1& v, const T2& c) {
return v * c;
}

template<typename T1, typename T2, int R1,int C1>
inline
Expand Down

0 comments on commit 2b0107e

Please sign in to comment.