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&eacute; approximant to the exponential.
60    *
61    *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
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&eacute; approximant to the exponential.
78    *
79    *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
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&eacute; approximant to the exponential.
96    *
97    *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
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&eacute; approximant to the exponential.
117    *
118    *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
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&eacute; approximant to the exponential.
140    *
141    *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
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&eacute; 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&eacute;
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&eacute; 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