1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2009, 2010, 2013 Jitse Niesen <jitse@maths.leeds.ac.uk> 5 // Copyright (C) 2011, 2013 Chen-Pang He <jdh8@ms63.hinet.net> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 // 11 // This file was edited for to the Stan math library to create 12 // the matrix exponential function (matrix_exp), 2016. 13 14 #ifndef STAN_MATH_PRIM_FUN_MATRIXEXPONENTIAL_H 15 #define STAN_MATH_PRIM_FUN_MATRIXEXPONENTIAL_H 16 17 #include <stan/math/prim/fun/value_of_rec.hpp> 18 #include <cmath> 19 20 namespace Eigen { 21 22 template <typename RealScalar> 23 struct MatrixExponentialScalingOp 24 { 25 /** \brief Constructor. 26 * 27 * \param[in] squarings The integer \f$ s \f$ in this document. 28 */ MatrixExponentialScalingOpMatrixExponentialScalingOp29 MatrixExponentialScalingOp(int squarings) : m_squarings(squarings) { } 30 31 32 /** \brief Scale a matrix coefficient. 33 * 34 * \param[in,out] x The scalar to be scaled, becoming \f$ 2^{-s} x \f$. 35 */ operatorMatrixExponentialScalingOp36 inline const RealScalar operator() (const RealScalar& x) const 37 { 38 using std::ldexp; 39 return ldexp(x, -m_squarings); 40 } 41 42 using ComplexScalar = std::complex<RealScalar>; 43 44 /** \brief Scale a matrix coefficient. 45 * 46 * \param[in,out] x The scalar to be scaled, becoming \f$ 2^{-s} x \f$. 47 */ operatorMatrixExponentialScalingOp48 inline const ComplexScalar operator() (const ComplexScalar& x) const 49 { 50 using std::ldexp; 51 return ComplexScalar(ldexp(x.real(), -m_squarings), ldexp(x.imag(), -m_squarings)); 52 } 53 54 private: 55 int m_squarings; 56 }; 57 58 59 /** \brief Compute the (3,3)-Padé approximant to the exponential. 60 * 61 * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé 62 * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. 63 */ 64 template <typename MatrixType> matrix_exp_pade3(const MatrixType & A,MatrixType & U,MatrixType & V)65 void matrix_exp_pade3(const MatrixType &A, MatrixType &U, MatrixType &V) 66 { 67 using Eigen::internal::traits; 68 using RealScalar = 69 typename Eigen::NumTraits<typename traits<MatrixType>::Scalar>::Real; 70 const RealScalar b[] = {120.L, 60.L, 12.L, 1.L}; 71 const MatrixType A2 = A * A; 72 const MatrixType tmp = b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols()); 73 U.noalias() = A * tmp; 74 V = b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols()); 75 } 76 77 /** \brief Compute the (5,5)-Padé approximant to the exponential. 78 * 79 * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé 80 * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. 81 */ 82 template <typename MatrixType> matrix_exp_pade5(const MatrixType & A,MatrixType & U,MatrixType & V)83 void matrix_exp_pade5(const MatrixType &A, MatrixType &U, MatrixType &V) 84 { 85 using RealScalar = typename Eigen::NumTraits< 86 typename Eigen::internal::traits<MatrixType>::Scalar>::Real; 87 const RealScalar b[] = {30240.L, 15120.L, 3360.L, 420.L, 30.L, 1.L}; 88 const MatrixType A2 = A * A; 89 const MatrixType A4 = A2 * A2; 90 const MatrixType tmp = b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols()); 91 U.noalias() = A * tmp; 92 V = b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols()); 93 } 94 95 /** \brief Compute the (7,7)-Padé approximant to the exponential. 96 * 97 * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé 98 * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. 99 */ 100 template <typename MatrixType> matrix_exp_pade7(const MatrixType & A,MatrixType & U,MatrixType & V)101 void matrix_exp_pade7(const MatrixType &A, MatrixType &U, MatrixType &V) 102 { 103 using Eigen::internal::traits; 104 using RealScalar = 105 typename Eigen::NumTraits<typename traits<MatrixType>::Scalar>::Real; 106 const RealScalar b[] = {17297280.L, 8648640.L, 1995840.L, 277200.L, 25200.L, 1512.L, 56.L, 1.L}; 107 const MatrixType A2 = A * A; 108 const MatrixType A4 = A2 * A2; 109 const MatrixType A6 = A4 * A2; 110 const MatrixType tmp = b[7] * A6 + b[5] * A4 + b[3] * A2 111 + b[1] * MatrixType::Identity(A.rows(), A.cols()); 112 U.noalias() = A * tmp; 113 V = b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols()); 114 } 115 116 /** \brief Compute the (9,9)-Padé approximant to the exponential. 117 * 118 * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé 119 * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. 120 */ 121 template <typename MatrixType> matrix_exp_pade9(const MatrixType & A,MatrixType & U,MatrixType & V)122 void matrix_exp_pade9(const MatrixType &A, MatrixType &U, MatrixType &V) 123 { 124 using Eigen::internal::traits; 125 using RealScalar = 126 typename Eigen::NumTraits<typename traits<MatrixType>::Scalar>::Real; 127 const RealScalar b[] = {17643225600.L, 8821612800.L, 2075673600.L, 302702400.L, 30270240.L, 128 2162160.L, 110880.L, 3960.L, 90.L, 1.L}; 129 const MatrixType A2 = A * A; 130 const MatrixType A4 = A2 * A2; 131 const MatrixType A6 = A4 * A2; 132 const MatrixType A8 = A6 * A2; 133 const MatrixType tmp = b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2 134 + b[1] * MatrixType::Identity(A.rows(), A.cols()); 135 U.noalias() = A * tmp; 136 V = b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols()); 137 } 138 139 /** \brief Compute the (13,13)-Padé approximant to the exponential. 140 * 141 * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé 142 * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. 143 */ 144 template <typename MatrixType> matrix_exp_pade13(const MatrixType & A,MatrixType & U,MatrixType & V)145 void matrix_exp_pade13(const MatrixType &A, MatrixType &U, MatrixType &V) 146 { 147 using Eigen::internal::traits; 148 using RealScalar = 149 typename Eigen::NumTraits<typename traits<MatrixType>::Scalar>::Real; 150 const RealScalar b[] = {64764752532480000.L, 32382376266240000.L, 7771770303897600.L, 151 1187353796428800.L, 129060195264000.L, 10559470521600.L, 670442572800.L, 152 33522128640.L, 1323241920.L, 40840800.L, 960960.L, 16380.L, 182.L, 1.L}; 153 const MatrixType A2 = A * A; 154 const MatrixType A4 = A2 * A2; 155 const MatrixType A6 = A4 * A2; 156 V = b[13] * A6 + b[11] * A4 + b[9] * A2; // used for temporary storage 157 MatrixType tmp = A6 * V; 158 tmp += b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols()); 159 U.noalias() = A * tmp; 160 tmp = b[12] * A6 + b[10] * A4 + b[8] * A2; 161 V.noalias() = A6 * tmp; 162 V += b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols()); 163 } 164 165 template <typename MatrixType, typename RealScalar = typename Eigen:: NumTraits<typename Eigen::internal::traits<MatrixType>::Scalar>::Real> 166 struct matrix_exp_computeUV 167 { 168 /** \brief Compute Padé approximant to the exponential. 169 * 170 * Computes \c U, \c V and \c squarings such that \f$ (V+U)(V-U)^{-1} \f$ is a Padé 171 * approximant of \f$ \exp(2^{-\mbox{squarings}}M) \f$ around \f$ M = 0 \f$, where \f$ M \f$ 172 * denotes the matrix \c arg. The degree of the Padé approximant and the value of squarings 173 * are chosen such that the approximation error is no more than the round-off error. 174 * 175 * <p> Edit for Stan: template ComputeUV::run so that it may used on 176 * autodiff variables (var and fvar). This required adding the scalar_type 177 * argument, which tells the function the type of elements used in the 178 * matrix. 179 */ 180 static void run(const MatrixType& arg, MatrixType& U, MatrixType& V, int& squarings); 181 }; 182 183 template <typename MatrixType> 184 struct matrix_exp_computeUV<MatrixType> 185 { 186 template <typename T> 187 static void run(const MatrixType& arg, MatrixType& U, MatrixType& V, int& squarings, 188 T scalar_type) 189 { 190 using std::frexp; 191 using std::pow; 192 using stan::math::value_of_rec; 193 const T l1norm = arg.cwiseAbs().colwise().sum().maxCoeff(); 194 squarings = 0; 195 if (l1norm < 1.495585217958292e-002) { 196 matrix_exp_pade3(arg, U, V); 197 } else if (l1norm < 2.539398330063230e-001) { 198 matrix_exp_pade5(arg, U, V); 199 } else if (l1norm < 9.504178996162932e-001) { 200 matrix_exp_pade7(arg, U, V); 201 } else if (l1norm < 2.097847961257068e+000) { 202 matrix_exp_pade9(arg, U, V); 203 } else { 204 const double maxnorm = 5.371920351148152; 205 frexp(value_of_rec(l1norm) / value_of_rec(maxnorm), &squarings); 206 if (squarings < 0) { 207 squarings = 0; 208 } 209 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<T>(squarings)); 210 matrix_exp_pade13(A, U, V); 211 } 212 } 213 }; 214 215 } 216 217 #endif 218