Stan  1.0
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
matrix.hpp
Go to the documentation of this file.
1 #ifndef __STAN__MATH__MATRIX_HPP__
2 #define __STAN__MATH__MATRIX_HPP__
3 
4 #include <stdexcept>
5 #include <ostream>
6 #include <vector>
7 
8 #include <boost/math/tools/promotion.hpp>
9 
10 #define EIGEN_DENSEBASE_PLUGIN "stan/math/EigenDenseBaseAddons.hpp"
11 #include <Eigen/Dense>
12 
13 #include <stan/math/boost_error_handling.hpp>
14 
15 namespace stan {
16 
17  namespace math {
18 
24  template<typename T>
25  struct EigenType {
26 
29  typedef T scalar;
30 
34  typedef Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> matrix;
35 
39  typedef Eigen::Matrix<T,Eigen::Dynamic,1> vector;
40 
44  typedef Eigen::Matrix<T,1,Eigen::Dynamic> row_vector;
45 
49  typedef Eigen::DiagonalMatrix<T,Eigen::Dynamic> diagonal_matrix;
50 
54  typedef Eigen::Array<T,1,Eigen::Dynamic> array_vector;
55 
59  typedef Eigen::Array<T,Eigen::Dynamic,1> array_row_vector;
60 
64  typedef Eigen::Array<T,Eigen::Dynamic,Eigen::Dynamic> array_matrix;
65  };
66 
70  typedef EigenType<double>::matrix matrix_d;
71 
75  typedef EigenType<double>::vector vector_d;
76 
80  typedef EigenType<double>::row_vector row_vector_d;
81 
82  namespace {
83 
84  template <typename T>
85  void resize(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& x,
86  const std::vector<size_t>& dims,
87  size_t pos) {
88  x.resize(dims[pos],dims[pos+1]);
89  }
90 
91  template <typename T>
92  void resize(Eigen::Matrix<T,Eigen::Dynamic,1>& x,
93  const std::vector<size_t>& dims,
94  size_t pos) {
95  x.resize(dims[pos]);
96  }
97 
98  template <typename T>
99  void resize(Eigen::Matrix<T,1,Eigen::Dynamic>& x,
100  const std::vector<size_t>& dims,
101  size_t pos) {
102  x.resize(dims[pos]);
103  }
104 
105 
106  void resize(double x,
107  const std::vector<size_t>& dims,
108  size_t pos) {
109  // no-op
110  }
111 
112  template <typename T>
113  void resize(std::vector<T>& x,
114  const std::vector<size_t>& dims,
115  size_t pos) {
116  x.resize(dims[pos]);
117  ++pos;
118  if (pos >= dims.size()) return; // skips lowest loop to scalar
119  for (size_t i = 0; i < x.size(); ++i)
120  resize(x[i],dims,pos);
121  }
122 
123  }
124 
134  template <typename T>
135  inline void resize(T& x, std::vector<size_t> dims) {
136  resize(x,dims,0U);
137  }
138 
139  // polymorphic gets with bounds checking
140 
141 
142  namespace {
143 
144  void raise_range_error(size_t max,
145  size_t i,
146  const char* msg,
147  size_t idx) {
148  std::stringstream s;
149  s << "INDEX OPERATOR [] OUT OF BOUNDS"
150  << "; index=" << i
151  << "; lower bound=1"
152  << "; upper bound=" << max
153  << "; index position=" << idx
154  << "; " << msg
155  << std::endl;
156  throw std::out_of_range(s.str());
157  }
158 
159  inline
160  void check_range(size_t max,
161  size_t i,
162  const char* msg,
163  size_t idx) {
164 #ifndef NDEBUG
165  if (i < 1 || i > max)
166  raise_range_error(max,i,msg,idx);
167 #endif
168  }
169 
170  }
171 
186  template <typename T>
187  inline
188  T& get_base1(std::vector<T>& x,
189  size_t i,
190  const char* error_msg,
191  size_t idx) {
192  check_range(x.size(),i,error_msg,idx);
193  return x[i - 1];
194  }
195 
211  template <typename T>
212  inline
213  T& get_base1(std::vector<std::vector<T> >& x,
214  size_t i1,
215  size_t i2,
216  const char* error_msg,
217  size_t idx) {
218  check_range(x.size(),i1,error_msg,idx);
219  return get_base1(x[i1 - 1],i2,error_msg,idx+1);
220  }
221 
238  template <typename T>
239  inline
240  T& get_base1(std::vector<std::vector<std::vector<T> > >& x,
241  size_t i1,
242  size_t i2,
243  size_t i3,
244  const char* error_msg,
245  size_t idx) {
246  check_range(x.size(),i1,error_msg,idx);
247  return get_base1(x[i1 - 1],i2,i3,error_msg,idx+1);
248  }
249 
267  template <typename T>
268  inline
269  T& get_base1(std::vector<std::vector<std::vector<std::vector<T> > > >& x,
270  size_t i1,
271  size_t i2,
272  size_t i3,
273  size_t i4,
274  const char* error_msg,
275  size_t idx) {
276  check_range(x.size(),i1,error_msg,idx);
277  return get_base1(x[i1 - 1],i2,i3,i4,error_msg,idx+1);
278  }
279 
298  template <typename T>
299  inline
300  T& get_base1(std::vector<std::vector<std::vector<std::vector<std::vector<T> > > > >& x,
301  size_t i1,
302  size_t i2,
303  size_t i3,
304  size_t i4,
305  size_t i5,
306  const char* error_msg,
307  size_t idx) {
308  check_range(x.size(),i1,error_msg,idx);
309  return get_base1(x[i1 - 1],i2,i3,i4,i5,error_msg,idx+1);
310  }
311 
331  template <typename T>
332  inline
333  T& get_base1(std::vector<std::vector<std::vector<std::vector<std::vector<std::vector<T> > > > > >& x,
334  size_t i1,
335  size_t i2,
336  size_t i3,
337  size_t i4,
338  size_t i5,
339  size_t i6,
340  const char* error_msg,
341  size_t idx) {
342  check_range(x.size(),i1,error_msg,idx);
343  return get_base1(x[i1 - 1],i2,i3,i4,i5,i6,error_msg,idx+1);
344  }
345 
346 
367  template <typename T>
368  inline
369  T& get_base1(std::vector<std::vector<std::vector<std::vector<std::vector<std::vector<std::vector<T> > > > > > >& x,
370  size_t i1,
371  size_t i2,
372  size_t i3,
373  size_t i4,
374  size_t i5,
375  size_t i6,
376  size_t i7,
377  const char* error_msg,
378  size_t idx) {
379  check_range(x.size(),i1,error_msg,idx);
380  return get_base1(x[i1 - 1],i2,i3,i4,i5,i6,i7,error_msg,idx+1);
381  }
382 
383 
405  template <typename T>
406  inline
407  T& get_base1(std::vector<std::vector<std::vector<std::vector<std::vector<std::vector<std::vector<std::vector<T> > > > > > > >& x,
408  size_t i1,
409  size_t i2,
410  size_t i3,
411  size_t i4,
412  size_t i5,
413  size_t i6,
414  size_t i7,
415  size_t i8,
416  const char* error_msg,
417  size_t idx) {
418  check_range(x.size(),i1,error_msg,idx);
419  return get_base1(x[i1 - 1],i2,i3,i4,i5,i6,i7,i8,error_msg,idx+1);
420  }
421 
422 
423 
443  template <typename T>
444  inline
445  Eigen::Matrix<T,1,Eigen::Dynamic>
446  get_base1(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& x,
447  size_t m,
448  const char* error_msg,
449  size_t idx) {
450  check_range(x.rows(),m,error_msg,idx);
451  return x.row(m - 1);
452  }
453 
470  template <typename T>
471  inline
472  T& get_base1(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& x,
473  size_t m,
474  size_t n,
475  const char* error_msg,
476  size_t idx) {
477  check_range(x.rows(),m,error_msg,idx);
478  check_range(x.cols(),n,error_msg,idx + 1);
479  return x(m - 1, n - 1);
480  }
481 
496  template <typename T>
497  inline
498  T& get_base1(Eigen::Matrix<T,Eigen::Dynamic,1>& x,
499  size_t m,
500  const char* error_msg,
501  size_t idx) {
502  check_range(x.size(),m,error_msg,idx);
503  return x(m - 1);
504 
505  }
506 
521  template <typename T>
522  inline
523  T& get_base1(Eigen::Matrix<T,1,Eigen::Dynamic>& x,
524  size_t n,
525  const char* error_msg,
526  size_t idx) {
527  check_range(x.size(),n,error_msg,idx);
528  return x(n - 1);
529  }
530 
531 
532  // int returns
533 
540  inline size_t rows(const vector_d& v) {
541  return v.size();
542  }
549  inline size_t rows(const row_vector_d& rv) {
550  return 1;
551  }
558  inline size_t rows(const matrix_d& m) {
559  return m.rows();
560  }
561 
568  inline size_t cols(const vector_d& v) {
569  return 1;
570  }
577  inline size_t cols(const row_vector_d& rv) {
578  return rv.size();
579  }
585  inline size_t cols(const matrix_d& m) {
586  return m.cols();
587  }
588 
589  template <typename T, int R, int C>
590  void validate_column_index(const Eigen::Matrix<T,R,C>& m,
591  size_t j,
592  const char* msg) {
593  if (j > 0 && j <= static_cast<size_t>(m.cols())) return;
594  std::stringstream ss;
595  ss << "require 0 < column index <= number of columns in" << msg;
596  ss << " found cols()=" << m.cols()
597  << "; index j=" << j;
598  throw std::domain_error(ss.str());
599  }
600 
601  template <typename T, int R, int C>
602  void validate_row_index(const Eigen::Matrix<T,R,C>& m,
603  size_t i,
604  const char* msg) {
605  if (i > 0 && i <= static_cast<size_t>(m.rows())) return;
606  std::stringstream ss;
607  ss << "require 0 < row index <= number of rows in" << msg;
608  ss << " found rows()=" << m.rows()
609  << "; index i=" << i;
610  throw std::domain_error(ss.str());
611  }
612 
613  template <typename T, int R, int C>
614  void validate_square(const Eigen::Matrix<T,R,C>& x,
615  const char* msg) {
616  if (x.rows() == x.cols()) return;
617  std::stringstream ss;
618  ss << "error in call to " << msg
619  << "; require square matrix, but found"
620  << " rows=" << x.rows()
621  << "; cols=" << x.cols();
622  throw std::domain_error(ss.str());
623  }
624 
625  template <typename T1, int R1, int C1, typename T2, int R2, int C2>
626  inline void validate_matching_dims(const Eigen::Matrix<T1,R1,C1>& x1,
627  const Eigen::Matrix<T2,R2,C2>& x2,
628  const char* msg) {
629  if (x1.rows() == x2.rows()
630  && x1.cols() == x2.cols()) return;
631  std::stringstream ss;
632  ss << "error in call to " << msg
633  << "; require matching dimensions, but found"
634  << " arg1 rows=" << x1.rows() << " arg1 cols=" << x1.cols()
635  << " arg2 rows=" << x2.rows() << " arg2 cols=" << x2.cols();
636  throw std::domain_error(ss.str());
637  }
638 
639  template <typename T1, typename T2>
640  inline void validate_matching_sizes(const std::vector<T1>& x1,
641  const std::vector<T2>& x2,
642  const char* msg) {
643  if (x1.size() == x2.size()) return;
644  std::stringstream ss;
645  ss << "require matching sizes in " << msg
646  << " found first argument size=" << x1.size()
647  << "; second argument size=" << x2.size();
648  throw std::domain_error(ss.str());
649  }
650 
651  template <typename T1, int R1, int C1, typename T2, int R2, int C2>
652  inline void validate_matching_sizes(const Eigen::Matrix<T1,R1,C1>& x1,
653  const Eigen::Matrix<T2,R2,C2>& x2,
654  const char* msg) {
655  if (x1.size() == x2.size()) return;
656  std::stringstream ss;
657  ss << "error in call to " << msg
658  << "; require matching sizes, but found"
659  << " arg1 rows=" << x1.rows() << " arg1 cols=" << x1.cols()
660  << " arg1 size=" << (x1.rows() * x1.cols())
661  << " arg2 rows=" << x2.rows() << " arg2 cols=" << x2.cols()
662  << " arg2 size=" << (x2.rows() * x2.cols());
663  throw std::domain_error(ss.str());
664  }
665 
666  template <typename T1, int R1, int C1, typename T2, int R2, int C2>
667  inline void validate_multiplicable(const Eigen::Matrix<T1,R1,C1>& x1,
668  const Eigen::Matrix<T2,R2,C2>& x2,
669  const char* msg) {
670  if (x1.cols() == x2.rows()) return;
671  std::stringstream ss;
672  ss << "error in call to " << msg
673  << "; require cols of arg1 to match rows of arg2, but found "
674  << " arg1 rows=" << x1.rows() << " arg1 cols=" << x1.cols()
675  << " arg2 rows=" << x2.rows() << " arg2 cols=" << x2.cols();
676  throw std::domain_error(ss.str());
677  }
678 
679  template <typename T>
680  inline void validate_nonzero_size(const T& x, const char* msg) {
681  if (x.size() > 0) return;
682  std::stringstream ss;
683  ss << "require non-zero size for " << msg
684  << "found size=" << x.size();
685  throw std::domain_error(ss.str());
686  }
687 
688  template <typename T, int R, int C>
689  inline void validate_vector(const Eigen::Matrix<T,R,C>& x,
690  const char* msg) {
691  if (x.rows() == 1 || x.cols() == 1) return;
692  std::stringstream ss;
693  ss << "error in " << msg
694  << "; require vector, found "
695  << " rows=" << x.rows() << "cols=" << x.cols();
696  throw std::domain_error(ss.str());
697  }
698 
699 
700 
701 
702 
703 
704  // scalar returns
705 
712  double determinant(const matrix_d& m);
713 
714 
722  template <int R, int C>
723  inline double dot_self(const Eigen::Matrix<double, R, C>& v) {
724  validate_vector(v,"dot_self");
725  double sum = 0.0;
726  for (int i = 0; i < v.size(); ++i)
727  sum += v(i) * v(i);
728  return sum;
729  }
730 
736  template<typename T>
737  inline Eigen::Matrix<T,Eigen::Dynamic,1>
738  columns_dot_self(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& x) {
739  return x.colwise().squaredNorm();
740  }
741 
751  template<int R1,int C1,int R2, int C2>
752  inline double dot_product(const Eigen::Matrix<double, R1, C1>& v1,
753  const Eigen::Matrix<double, R2, C2>& v2) {
754  validate_vector(v1,"dot_product");
755  validate_vector(v2,"dot_product");
756  validate_matching_sizes(v1,v2,"dot_product");
757  double sum = 0.0;
758  for (int i = 0; i < v1.size(); ++i)
759  sum += v1[i] * v2[i];
760  return sum;
761  }
768  inline double dot_product(const double* v1, const double* v2,
769  size_t length) {
770  double result = 0;
771  for (size_t i = 0; i < length; i++)
772  result += v1[i] * v2[i];
773  return result;
774  }
781  inline double dot_product(const std::vector<double>& v1,
782  const std::vector<double>& v2) {
783  validate_matching_sizes(v1,v2,"dot_product");
784  return dot_product(&v1[0], &v2[0], v1.size());
785  }
786 
794  inline int min(const std::vector<int>& x) {
795  if (x.size() == 0)
796  throw std::domain_error("error: cannot take min of empty int vector");
797  int min = x[0];
798  for (size_t i = 1; i < x.size(); ++i)
799  if (x[i] < min)
800  min = x[i];
801  return min;
802  }
803 
811  template <typename T>
812  inline T min(const std::vector<T>& x) {
813  if (x.size() == 0)
814  return std::numeric_limits<T>::infinity();
815  T min = x[0];
816  for (size_t i = 1; i < x.size(); ++i)
817  if (x[i] < min)
818  min = x[i];
819  return min;
820  }
821 
828  template <typename T, int R, int C>
829  inline T min(const Eigen::Matrix<T,R,C>& m) {
830  if (m.size() == 0)
831  return std::numeric_limits<double>::infinity();
832  return m.minCoeff();
833  }
834 
835 
836 
845  inline int max(const std::vector<int>& x) {
846  if (x.size() == 0)
847  throw std::domain_error("error: cannot take max of empty int vector");
848  int max = x[0];
849  for (size_t i = 1; i < x.size(); ++i)
850  if (x[i] > max)
851  max = x[i];
852  return max;
853  }
854 
862  template <typename T>
863  inline T max(const std::vector<T>& x) {
864  if (x.size() == 0)
865  return -std::numeric_limits<T>::infinity();
866  T max = x[0];
867  for (size_t i = 1; i < x.size(); ++i)
868  if (x[i] > max)
869  max = x[i];
870  return max;
871  }
872 
879  template <typename T, int R, int C>
880  inline T max(const Eigen::Matrix<T,R,C>& m) {
881  if (m.size() == 0)
882  return -std::numeric_limits<double>::infinity();
883  return m.maxCoeff();
884  }
885 
894  template <typename T>
895  inline
896  typename boost::math::tools::promote_args<T>::type
897  mean(const std::vector<T>& v) {
898  validate_nonzero_size(v,"mean");
899  T sum(v[0]);
900  for (size_t i = 1; i < v.size(); ++i)
901  sum += v[i];
902  return sum / v.size();
903  }
904 
911  template <typename T, int R, int C>
912  inline
913  typename boost::math::tools::promote_args<T>::type
914  mean(const Eigen::Matrix<T,R,C>& m) {
915  validate_nonzero_size(m,"mean");
916  return m.mean();
917  }
918 
927  template <typename T>
928  inline
929  typename boost::math::tools::promote_args<T>::type
930  variance(const std::vector<T>& v) {
931  validate_nonzero_size(v,"variance");
932  if (v.size() == 1)
933  return 0.0;
934  T v_mean(mean(v));
935  T sum_sq_diff(0);
936  for (size_t i = 0; i < v.size(); ++i) {
937  T diff = v[i] - v_mean;
938  sum_sq_diff += diff * diff;
939  }
940  return sum_sq_diff / (v.size() - 1);
941  }
942 
949  template <typename T, int R, int C>
950  inline
951  typename boost::math::tools::promote_args<T>::type
952  variance(const Eigen::Matrix<T,R,C>& m) {
953  validate_nonzero_size(m,"variance");
954  if (m.size() == 1)
955  return 0.0;
956  typename boost::math::tools::promote_args<T>::type
957  mn(mean(m));
958  typename boost::math::tools::promote_args<T>::type
959  sum_sq_diff(0);
960  for (int i = 0; i < m.size(); ++i) {
961  typename boost::math::tools::promote_args<T>::type
962  diff = m(i) - mn;
963  sum_sq_diff += diff * diff;
964  }
965  return sum_sq_diff / (m.size() - 1);
966  }
967 
974  template <typename T>
975  inline
976  typename boost::math::tools::promote_args<T>::type
977  sd(const std::vector<T>& v) {
978  validate_nonzero_size(v,"sd");
979  if (v.size() == 1) return 0.0;
980  return sqrt(variance(v));
981  }
982 
989  template <typename T, int R, int C>
990  inline
991  typename boost::math::tools::promote_args<T>::type
992  sd(const Eigen::Matrix<T,R,C>& m) {
993  // FIXME: redundant with test in variance; second line saves sqrt
994  validate_nonzero_size(m,"sd");
995  if (m.size() == 1) return 0.0;
996  return sqrt(variance(m));
997  }
998 
1007  template <typename T>
1008  inline T sum(const std::vector<T>& xs) {
1009  if (xs.size() == 0) return 0;
1010  T sum(xs[0]);
1011  for (size_t i = 1; i < xs.size(); ++i)
1012  sum += xs[i];
1013  return sum;
1014  }
1015 
1022  template <typename T, int R, int C>
1023  inline double sum(const Eigen::Matrix<T,R,C>& v) {
1024  return v.sum();
1025  }
1026 
1033  template <typename T>
1034  inline T prod(const std::vector<T>& v) {
1035  if (v.size() == 0) return 1;
1036  T product = v[0];
1037  for (size_t i = 1; i < v.size(); ++i)
1038  product *= v[i];
1039  return product;
1040  }
1041 
1048  template <typename T, int R, int C>
1049  inline T prod(const Eigen::Matrix<T,R,C>& v) {
1050  if (v.size() == 0) return 1.0;
1051  return v.prod();
1052  }
1053 
1063  template <typename T>
1064  inline T trace(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m) {
1065  return m.trace();
1066  }
1067 
1074  template<typename T, int Rows, int Cols>
1075  inline Eigen::Matrix<T,Rows,Cols> log(const Eigen::Matrix<T,Rows,Cols>& m) {
1076  return m.array().log().matrix();
1077  }
1078 
1085  template<typename T, int Rows, int Cols>
1086  inline Eigen::Matrix<T,Rows,Cols> exp(const Eigen::Matrix<T,Rows,Cols>& m) {
1087  return m.array().exp().matrix();
1088  }
1089 
1090 
1091  // vector and matrix returns
1092 
1101  vector_d add(const vector_d& v1, const vector_d& v2);
1110  row_vector_d add(const row_vector_d& rv1,
1111  const row_vector_d& rv2);
1120  matrix_d add(const matrix_d& m1, const matrix_d& m2);
1121 
1122 
1129  template<int Rows, int Cols>
1130  inline Eigen::Matrix<double,Rows,Cols> add(const Eigen::Matrix<double,Rows,Cols>& m, const double& c) {
1131  return (m.array() + c).matrix();
1132  }
1139  template<int Rows, int Cols>
1140  inline Eigen::Matrix<double,Rows,Cols> add(const double& c, const Eigen::Matrix<double,Rows,Cols>& m) {
1141  return (c + m.array()).matrix();
1142  }
1143 
1152  vector_d subtract(const vector_d& v1, const vector_d& v2);
1161  row_vector_d subtract(const row_vector_d& rv1, const row_vector_d& rv2);
1170  matrix_d subtract(const matrix_d& m1, const matrix_d& m2);
1171 
1178  template<int Rows, int Cols>
1179  inline Eigen::Matrix<double,Rows,Cols> subtract(const Eigen::Matrix<double,Rows,Cols>& m,
1180  double c) {
1181  return (m.array() - c).matrix();
1182  }
1189  template<int Rows, int Cols>
1190  inline Eigen::Matrix<double,Rows,Cols> subtract(double c, const Eigen::Matrix<double,Rows,Cols>& m) {
1191  return (c - m.array()).matrix();
1192  }
1193 
1200  vector_d minus(const vector_d& v);
1207  row_vector_d minus(const row_vector_d& rv);
1214  matrix_d minus(const matrix_d& m);
1215 
1223  vector_d divide(const vector_d& v, double c);
1231  row_vector_d divide(const row_vector_d& rv, double c);
1239  matrix_d divide(const matrix_d& m, double c);
1240 
1247  vector_d elt_multiply(const vector_d& v1, const vector_d& v2);
1254  row_vector_d elt_multiply(const row_vector_d& v1, const row_vector_d& v2);
1261  matrix_d elt_multiply(const matrix_d& m1, const matrix_d& m2);
1262 
1263 
1270  vector_d elt_divide(const vector_d& v1, const vector_d& v2);
1277  row_vector_d elt_divide(const row_vector_d& v1, const row_vector_d& v2);
1284  matrix_d elt_divide(const matrix_d& m1, const matrix_d& m2);
1285 
1293  vector_d multiply(const vector_d& v, double c);
1301  row_vector_d multiply(const row_vector_d& rv, double c);
1309  matrix_d multiply(const matrix_d& m, double c);
1320  template<int R1,int C1,int R2,int C2>
1321  inline Eigen::Matrix<double,R1,C2> multiply(const Eigen::Matrix<double,R1,C1>& m1,
1322  const Eigen::Matrix<double,R2,C2>& m2) {
1323 
1324  if (m1.cols() != m2.rows())
1325  throw std::domain_error("m1.cols() != m2.rows()");
1326  return m1*m2;
1327  }
1337  template<int C1,int R2>
1338  inline double multiply(const Eigen::Matrix<double,1,C1>& rv,
1339  const Eigen::Matrix<double,R2,1>& v) {
1340  stan::math::validate_matching_sizes(rv,v,"multiply");
1341  if (rv.size() != v.size())
1342  throw std::domain_error("rv.size() != v.size()");
1343  return rv.dot(v);
1344  }
1351  vector_d multiply(double c, const vector_d& v);
1358  row_vector_d multiply(double c, const row_vector_d& rv);
1365  matrix_d multiply(double c, const matrix_d& m);
1366 
1375  inline matrix_d
1376  multiply_lower_tri_self_transpose(const matrix_d& L) {
1377  stan::math::validate_square(L,"multiply_lower_tri_self_transpose");
1378  if (L.rows() == 0)
1379  return matrix_d(0,0);
1380  if (L.rows() == 1) {
1381  matrix_d result(1,1);
1382  result(0,0) = L(0,0) * L(0,0);
1383  return result;
1384  }
1385  // FIXME: write custom following agrad/matrix because can't get L_tri into
1386  // multiplication as no template support for tri * tri
1387  matrix_d L_tri = L.transpose().triangularView<Eigen::Upper>();
1388  return L.triangularView<Eigen::Lower>() * L_tri;
1389  }
1390 
1397  inline matrix_d
1398  tcrossprod(const matrix_d& M) {
1399  if(M.rows() == 0)
1400  return matrix_d(0,0);
1401  if(M.rows() == 1) {
1402  return M * M.transpose();
1403  }
1404  matrix_d result(M.rows(),M.rows());
1405  return result.setZero().selfadjointView<Eigen::Upper>().rankUpdate(M);
1406  }
1407 
1414  inline matrix_d
1415  crossprod(const matrix_d& M) {
1416  return tcrossprod(M.transpose());
1417  }
1418 
1430  row_vector_d row(const matrix_d& m, size_t i);
1431 
1443  vector_d col(const matrix_d& m, size_t j);
1444 
1451  vector_d diagonal(const matrix_d& m);
1452 
1459  matrix_d diag_matrix(const vector_d& v);
1460 
1467  row_vector_d transpose(const vector_d& v);
1474  vector_d transpose(const row_vector_d& rv);
1480  matrix_d transpose(const matrix_d& m);
1481 
1487  matrix_d inverse(const matrix_d& m);
1488 
1494  vector_d softmax(const vector_d& y);
1495 
1496 
1497  template<int R1,int C1,int R2,int C2>
1498  inline Eigen::Matrix<double,R1,C2> mdivide_left_tri_low(const Eigen::Matrix<double,R1,C1> &A,
1499  const Eigen::Matrix<double,R2,C2> &b) {
1500  if (A.cols() != A.rows())
1501  throw std::domain_error("A is not square");
1502  if (A.cols() != b.rows())
1503  throw std::domain_error("A.cols() != b.rows()");
1504  return A.template triangularView<Eigen::Lower>().solve(b);
1505  }
1506 
1507  inline Eigen::MatrixXd mdivide_left_tri_low(const Eigen::MatrixXd &A) {
1508  if (A.cols() != A.rows())
1509  throw std::domain_error("A is not square");
1510  int n = A.rows();
1511  Eigen::MatrixXd b;
1512  b.setIdentity(n,n);
1513  A.triangularView<Eigen::Lower>().solveInPlace(b);
1514  return b;
1515  }
1516 
1526  template<int TriView,int R1,int C1,int R2,int C2>
1527  inline Eigen::Matrix<double,R1,C2> mdivide_left_tri(const Eigen::Matrix<double,R1,C1> &A,
1528  const Eigen::Matrix<double,R2,C2> &b) {
1529  if (A.cols() != A.rows())
1530  throw std::domain_error("A is not square");
1531  if (A.cols() != b.rows())
1532  throw std::domain_error("A.cols() != b.rows()");
1533  return A.template triangularView<TriView>().solve(b);
1534  }
1542  template<int TriView>
1543  inline Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> mdivide_left_tri(const Eigen::MatrixXd &A) {
1544  if (A.cols() != A.rows())
1545  throw std::domain_error("A is not square");
1546  int n = A.rows();
1547  Eigen::MatrixXd b;
1548  b.setIdentity(n,n);
1549  A.triangularView<TriView>().solveInPlace(b);
1550  return b;
1551  }
1561  template<int TriView,int R1,int C1,int R2,int C2>
1562  inline Eigen::Matrix<double,R1,C2> mdivide_right_tri(const Eigen::Matrix<double,R1,C1> &b,
1563  const Eigen::Matrix<double,R2,C2> &A) {
1564  if (A.cols() != A.rows())
1565  throw std::domain_error("A is not square");
1566  if (A.rows() != b.cols())
1567  throw std::domain_error("A.rows() != b.cols()");
1568  return A.template triangularView<TriView>().transpose().solve(b.transpose()).transpose();
1569  }
1579  template<int R1,int C1,int R2,int C2>
1580  inline Eigen::Matrix<double,R1,C2> mdivide_right_tri_low(const Eigen::Matrix<double,R1,C1> &b,
1581  const Eigen::Matrix<double,R2,C2> &A) {
1582  return mdivide_right_tri<Eigen::Lower>(b,A);
1583  }
1584 
1585 
1594  template<int R1,int C1,int R2,int C2>
1595  inline Eigen::Matrix<double,R1,C2> mdivide_left(const Eigen::Matrix<double,R1,C1> &A,
1596  const Eigen::Matrix<double,R2,C2> &b) {
1597  if (A.cols() != A.rows())
1598  throw std::domain_error("A is not square");
1599  if (A.cols() != b.rows())
1600  throw std::domain_error("A.cols() != b.rows()");
1601  return A.lu().solve(b);
1602  }
1611  template<int R1,int C1,int R2,int C2>
1612  inline Eigen::Matrix<double,R1,C2> mdivide_right(const Eigen::Matrix<double,R1,C1> &b,
1613  const Eigen::Matrix<double,R2,C2> &A) {
1614  if (A.cols() != A.rows())
1615  throw std::domain_error("A is not square");
1616  if (A.rows() != b.cols())
1617  throw std::domain_error("A.rows() != b.cols()");
1618  return A.transpose().lu().solve(b.transpose()).transpose();
1619  }
1620 
1628  vector_d eigenvalues(const matrix_d& m);
1629 
1637  matrix_d eigenvectors(const matrix_d& m);
1655  void eigen_decompose(const matrix_d& m,
1656  vector_d& eigenvalues,
1657  matrix_d& eigenvectors);
1658 
1668  vector_d eigenvalues_sym(const matrix_d& m);
1678  matrix_d eigenvectors_sym(const matrix_d& m);
1691  void eigen_decompose_sym(const matrix_d& m,
1692  vector_d& eigenvalues,
1693  matrix_d& eigenvectors);
1694 
1695 
1706  matrix_d cholesky_decompose(const matrix_d& m);
1707 
1716  vector_d singular_values(const matrix_d& m);
1717 
1739  void svd(const matrix_d& m, matrix_d& u, matrix_d& v, vector_d& s);
1740 
1741  template <typename T>
1742  void stan_print(std::ostream* o, const T& x) {
1743  *o << x;
1744  }
1745 
1746  template <typename T>
1747  void stan_print(std::ostream* o, const std::vector<T>& x) {
1748  *o << '[';
1749  for (int i = 0; i < x.size(); ++i) {
1750  if (i > 0) *o << ',';
1751  stan_print(o,x[i]);
1752  }
1753  *o << ']';
1754  }
1755 
1756  template <typename T>
1757  void stan_print(std::ostream* o, const Eigen::Matrix<T,Eigen::Dynamic,1>& x) {
1758  *o << '[';
1759  for (int i = 0; i < x.size(); ++i) {
1760  if (i > 0) *o << ',';
1761  stan_print(o,x(i));
1762  }
1763  *o << ']';
1764  }
1765 
1766  template <typename T>
1767  void stan_print(std::ostream* o, const Eigen::Matrix<T,1,Eigen::Dynamic>& x) {
1768  *o << '[';
1769  for (int i = 0; i < x.size(); ++i) {
1770  if (i > 0) *o << ',';
1771  stan_print(o,x(i));
1772  }
1773  *o << ']';
1774  }
1775 
1776  template <typename T>
1777  void stan_print(std::ostream* o, const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& x) {
1778  *o << '[';
1779  for (int i = 0; i < x.rows(); ++i) {
1780  if (i > 0) *o << ',';
1781  stan_print(o,x.row(i));
1782  }
1783  *o << ']';
1784  }
1785 
1786 
1787  }
1788 
1789 }
1790 
1791 #endif
1792 
     [ Stan Home Page ] © 2011–2012, Stan Development Team.