• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009, 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
5 // Copyright (C) 2011 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 #ifndef EIGEN_MATRIX_EXPONENTIAL
12 #define EIGEN_MATRIX_EXPONENTIAL
13 
14 #include "StemFunction.h"
15 
16 namespace Eigen {
17 
18 /** \ingroup MatrixFunctions_Module
19   * \brief Class for computing the matrix exponential.
20   * \tparam MatrixType type of the argument of the exponential,
21   * expected to be an instantiation of the Matrix class template.
22   */
23 template <typename MatrixType>
24 class MatrixExponential {
25 
26   public:
27 
28     /** \brief Constructor.
29       *
30       * The class stores a reference to \p M, so it should not be
31       * changed (or destroyed) before compute() is called.
32       *
33       * \param[in] M  matrix whose exponential is to be computed.
34       */
35     MatrixExponential(const MatrixType &M);
36 
37     /** \brief Computes the matrix exponential.
38       *
39       * \param[out] result  the matrix exponential of \p M in the constructor.
40       */
41     template <typename ResultType>
42     void compute(ResultType &result);
43 
44   private:
45 
46     // Prevent copying
47     MatrixExponential(const MatrixExponential&);
48     MatrixExponential& operator=(const MatrixExponential&);
49 
50     /** \brief Compute the (3,3)-Pad&eacute; approximant to the exponential.
51      *
52      *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
53      *  approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
54      *
55      *  \param[in] A   Argument of matrix exponential
56      */
57     void pade3(const MatrixType &A);
58 
59     /** \brief Compute the (5,5)-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      *  \param[in] A   Argument of matrix exponential
65      */
66     void pade5(const MatrixType &A);
67 
68     /** \brief Compute the (7,7)-Pad&eacute; approximant to the exponential.
69      *
70      *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
71      *  approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
72      *
73      *  \param[in] A   Argument of matrix exponential
74      */
75     void pade7(const MatrixType &A);
76 
77     /** \brief Compute the (9,9)-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      *  \param[in] A   Argument of matrix exponential
83      */
84     void pade9(const MatrixType &A);
85 
86     /** \brief Compute the (13,13)-Pad&eacute; approximant to the exponential.
87      *
88      *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
89      *  approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
90      *
91      *  \param[in] A   Argument of matrix exponential
92      */
93     void pade13(const MatrixType &A);
94 
95     /** \brief Compute the (17,17)-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      *  This function activates only if your long double is double-double or quadruple.
101      *
102      *  \param[in] A   Argument of matrix exponential
103      */
104     void pade17(const MatrixType &A);
105 
106     /** \brief Compute Pad&eacute; approximant to the exponential.
107      *
108      * Computes \c m_U, \c m_V and \c m_squarings such that
109      * \f$ (V+U)(V-U)^{-1} \f$ is a Pad&eacute; of
110      * \f$ \exp(2^{-\mbox{squarings}}M) \f$ around \f$ M = 0 \f$. The
111      * degree of the Pad&eacute; approximant and the value of
112      * squarings are chosen such that the approximation error is no
113      * more than the round-off error.
114      *
115      * The argument of this function should correspond with the (real
116      * part of) the entries of \c m_M.  It is used to select the
117      * correct implementation using overloading.
118      */
119     void computeUV(double);
120 
121     /** \brief Compute Pad&eacute; approximant to the exponential.
122      *
123      *  \sa computeUV(double);
124      */
125     void computeUV(float);
126 
127     /** \brief Compute Pad&eacute; approximant to the exponential.
128      *
129      *  \sa computeUV(double);
130      */
131     void computeUV(long double);
132 
133     typedef typename internal::traits<MatrixType>::Scalar Scalar;
134     typedef typename NumTraits<Scalar>::Real RealScalar;
135     typedef typename std::complex<RealScalar> ComplexScalar;
136 
137     /** \brief Reference to matrix whose exponential is to be computed. */
138     typename internal::nested<MatrixType>::type m_M;
139 
140     /** \brief Odd-degree terms in numerator of Pad&eacute; approximant. */
141     MatrixType m_U;
142 
143     /** \brief Even-degree terms in numerator of Pad&eacute; approximant. */
144     MatrixType m_V;
145 
146     /** \brief Used for temporary storage. */
147     MatrixType m_tmp1;
148 
149     /** \brief Used for temporary storage. */
150     MatrixType m_tmp2;
151 
152     /** \brief Identity matrix of the same size as \c m_M. */
153     MatrixType m_Id;
154 
155     /** \brief Number of squarings required in the last step. */
156     int m_squarings;
157 
158     /** \brief L1 norm of m_M. */
159     RealScalar m_l1norm;
160 };
161 
162 template <typename MatrixType>
MatrixExponential(const MatrixType & M)163 MatrixExponential<MatrixType>::MatrixExponential(const MatrixType &M) :
164   m_M(M),
165   m_U(M.rows(),M.cols()),
166   m_V(M.rows(),M.cols()),
167   m_tmp1(M.rows(),M.cols()),
168   m_tmp2(M.rows(),M.cols()),
169   m_Id(MatrixType::Identity(M.rows(), M.cols())),
170   m_squarings(0),
171   m_l1norm(M.cwiseAbs().colwise().sum().maxCoeff())
172 {
173   /* empty body */
174 }
175 
176 template <typename MatrixType>
177 template <typename ResultType>
compute(ResultType & result)178 void MatrixExponential<MatrixType>::compute(ResultType &result)
179 {
180 #if LDBL_MANT_DIG > 112 // rarely happens
181   if(sizeof(RealScalar) > 14) {
182     result = m_M.matrixFunction(StdStemFunctions<ComplexScalar>::exp);
183     return;
184   }
185 #endif
186   computeUV(RealScalar());
187   m_tmp1 = m_U + m_V;   // numerator of Pade approximant
188   m_tmp2 = -m_U + m_V;  // denominator of Pade approximant
189   result = m_tmp2.partialPivLu().solve(m_tmp1);
190   for (int i=0; i<m_squarings; i++)
191     result *= result;   // undo scaling by repeated squaring
192 }
193 
194 template <typename MatrixType>
pade3(const MatrixType & A)195 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade3(const MatrixType &A)
196 {
197   const RealScalar b[] = {120., 60., 12., 1.};
198   m_tmp1.noalias() = A * A;
199   m_tmp2 = b[3]*m_tmp1 + b[1]*m_Id;
200   m_U.noalias() = A * m_tmp2;
201   m_V = b[2]*m_tmp1 + b[0]*m_Id;
202 }
203 
204 template <typename MatrixType>
pade5(const MatrixType & A)205 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade5(const MatrixType &A)
206 {
207   const RealScalar b[] = {30240., 15120., 3360., 420., 30., 1.};
208   MatrixType A2 = A * A;
209   m_tmp1.noalias() = A2 * A2;
210   m_tmp2 = b[5]*m_tmp1 + b[3]*A2 + b[1]*m_Id;
211   m_U.noalias() = A * m_tmp2;
212   m_V = b[4]*m_tmp1 + b[2]*A2 + b[0]*m_Id;
213 }
214 
215 template <typename MatrixType>
pade7(const MatrixType & A)216 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade7(const MatrixType &A)
217 {
218   const RealScalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.};
219   MatrixType A2 = A * A;
220   MatrixType A4 = A2 * A2;
221   m_tmp1.noalias() = A4 * A2;
222   m_tmp2 = b[7]*m_tmp1 + b[5]*A4 + b[3]*A2 + b[1]*m_Id;
223   m_U.noalias() = A * m_tmp2;
224   m_V = b[6]*m_tmp1 + b[4]*A4 + b[2]*A2 + b[0]*m_Id;
225 }
226 
227 template <typename MatrixType>
pade9(const MatrixType & A)228 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade9(const MatrixType &A)
229 {
230   const RealScalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240.,
231 		      2162160., 110880., 3960., 90., 1.};
232   MatrixType A2 = A * A;
233   MatrixType A4 = A2 * A2;
234   MatrixType A6 = A4 * A2;
235   m_tmp1.noalias() = A6 * A2;
236   m_tmp2 = b[9]*m_tmp1 + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*m_Id;
237   m_U.noalias() = A * m_tmp2;
238   m_V = b[8]*m_tmp1 + b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*m_Id;
239 }
240 
241 template <typename MatrixType>
pade13(const MatrixType & A)242 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade13(const MatrixType &A)
243 {
244   const RealScalar b[] = {64764752532480000., 32382376266240000., 7771770303897600.,
245 		      1187353796428800., 129060195264000., 10559470521600., 670442572800.,
246 		      33522128640., 1323241920., 40840800., 960960., 16380., 182., 1.};
247   MatrixType A2 = A * A;
248   MatrixType A4 = A2 * A2;
249   m_tmp1.noalias() = A4 * A2;
250   m_V = b[13]*m_tmp1 + b[11]*A4 + b[9]*A2; // used for temporary storage
251   m_tmp2.noalias() = m_tmp1 * m_V;
252   m_tmp2 += b[7]*m_tmp1 + b[5]*A4 + b[3]*A2 + b[1]*m_Id;
253   m_U.noalias() = A * m_tmp2;
254   m_tmp2 = b[12]*m_tmp1 + b[10]*A4 + b[8]*A2;
255   m_V.noalias() = m_tmp1 * m_tmp2;
256   m_V += b[6]*m_tmp1 + b[4]*A4 + b[2]*A2 + b[0]*m_Id;
257 }
258 
259 #if LDBL_MANT_DIG > 64
260 template <typename MatrixType>
pade17(const MatrixType & A)261 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade17(const MatrixType &A)
262 {
263   const RealScalar b[] = {830034394580628357120000.L, 415017197290314178560000.L,
264 		      100610229646136770560000.L, 15720348382208870400000.L,
265 		      1774878043152614400000.L, 153822763739893248000.L, 10608466464820224000.L,
266 		      595373117923584000.L, 27563570274240000.L, 1060137318240000.L,
267 		      33924394183680.L, 899510451840.L, 19554575040.L, 341863200.L, 4651200.L,
268 		      46512.L, 306.L, 1.L};
269   MatrixType A2 = A * A;
270   MatrixType A4 = A2 * A2;
271   MatrixType A6 = A4 * A2;
272   m_tmp1.noalias() = A4 * A4;
273   m_V = b[17]*m_tmp1 + b[15]*A6 + b[13]*A4 + b[11]*A2; // used for temporary storage
274   m_tmp2.noalias() = m_tmp1 * m_V;
275   m_tmp2 += b[9]*m_tmp1 + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*m_Id;
276   m_U.noalias() = A * m_tmp2;
277   m_tmp2 = b[16]*m_tmp1 + b[14]*A6 + b[12]*A4 + b[10]*A2;
278   m_V.noalias() = m_tmp1 * m_tmp2;
279   m_V += b[8]*m_tmp1 + b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*m_Id;
280 }
281 #endif
282 
283 template <typename MatrixType>
computeUV(float)284 void MatrixExponential<MatrixType>::computeUV(float)
285 {
286   using std::frexp;
287   using std::pow;
288   if (m_l1norm < 4.258730016922831e-001) {
289     pade3(m_M);
290   } else if (m_l1norm < 1.880152677804762e+000) {
291     pade5(m_M);
292   } else {
293     const float maxnorm = 3.925724783138660f;
294     frexp(m_l1norm / maxnorm, &m_squarings);
295     if (m_squarings < 0) m_squarings = 0;
296     MatrixType A = m_M / pow(Scalar(2), m_squarings);
297     pade7(A);
298   }
299 }
300 
301 template <typename MatrixType>
computeUV(double)302 void MatrixExponential<MatrixType>::computeUV(double)
303 {
304   using std::frexp;
305   using std::pow;
306   if (m_l1norm < 1.495585217958292e-002) {
307     pade3(m_M);
308   } else if (m_l1norm < 2.539398330063230e-001) {
309     pade5(m_M);
310   } else if (m_l1norm < 9.504178996162932e-001) {
311     pade7(m_M);
312   } else if (m_l1norm < 2.097847961257068e+000) {
313     pade9(m_M);
314   } else {
315     const double maxnorm = 5.371920351148152;
316     frexp(m_l1norm / maxnorm, &m_squarings);
317     if (m_squarings < 0) m_squarings = 0;
318     MatrixType A = m_M / pow(Scalar(2), m_squarings);
319     pade13(A);
320   }
321 }
322 
323 template <typename MatrixType>
computeUV(long double)324 void MatrixExponential<MatrixType>::computeUV(long double)
325 {
326   using std::frexp;
327   using std::pow;
328 #if   LDBL_MANT_DIG == 53   // double precision
329   computeUV(double());
330 #elif LDBL_MANT_DIG <= 64   // extended precision
331   if (m_l1norm < 4.1968497232266989671e-003L) {
332     pade3(m_M);
333   } else if (m_l1norm < 1.1848116734693823091e-001L) {
334     pade5(m_M);
335   } else if (m_l1norm < 5.5170388480686700274e-001L) {
336     pade7(m_M);
337   } else if (m_l1norm < 1.3759868875587845383e+000L) {
338     pade9(m_M);
339   } else {
340     const long double maxnorm = 4.0246098906697353063L;
341     frexp(m_l1norm / maxnorm, &m_squarings);
342     if (m_squarings < 0) m_squarings = 0;
343     MatrixType A = m_M / pow(Scalar(2), m_squarings);
344     pade13(A);
345   }
346 #elif LDBL_MANT_DIG <= 106  // double-double
347   if (m_l1norm < 3.2787892205607026992947488108213e-005L) {
348     pade3(m_M);
349   } else if (m_l1norm < 6.4467025060072760084130906076332e-003L) {
350     pade5(m_M);
351   } else if (m_l1norm < 6.8988028496595374751374122881143e-002L) {
352     pade7(m_M);
353   } else if (m_l1norm < 2.7339737518502231741495857201670e-001L) {
354     pade9(m_M);
355   } else if (m_l1norm < 1.3203382096514474905666448850278e+000L) {
356     pade13(m_M);
357   } else {
358     const long double maxnorm = 3.2579440895405400856599663723517L;
359     frexp(m_l1norm / maxnorm, &m_squarings);
360     if (m_squarings < 0) m_squarings = 0;
361     MatrixType A = m_M / pow(Scalar(2), m_squarings);
362     pade17(A);
363   }
364 #elif LDBL_MANT_DIG <= 112  // quadruple precison
365   if (m_l1norm < 1.639394610288918690547467954466970e-005L) {
366     pade3(m_M);
367   } else if (m_l1norm < 4.253237712165275566025884344433009e-003L) {
368     pade5(m_M);
369   } else if (m_l1norm < 5.125804063165764409885122032933142e-002L) {
370     pade7(m_M);
371   } else if (m_l1norm < 2.170000765161155195453205651889853e-001L) {
372     pade9(m_M);
373   } else if (m_l1norm < 1.125358383453143065081397882891878e+000L) {
374     pade13(m_M);
375   } else {
376     const long double maxnorm = 2.884233277829519311757165057717815L;
377     frexp(m_l1norm / maxnorm, &m_squarings);
378     if (m_squarings < 0) m_squarings = 0;
379     MatrixType A = m_M / pow(Scalar(2), m_squarings);
380     pade17(A);
381   }
382 #else
383   // this case should be handled in compute()
384   eigen_assert(false && "Bug in MatrixExponential");
385 #endif  // LDBL_MANT_DIG
386 }
387 
388 /** \ingroup MatrixFunctions_Module
389   *
390   * \brief Proxy for the matrix exponential of some matrix (expression).
391   *
392   * \tparam Derived  Type of the argument to the matrix exponential.
393   *
394   * This class holds the argument to the matrix exponential until it
395   * is assigned or evaluated for some other reason (so the argument
396   * should not be changed in the meantime). It is the return type of
397   * MatrixBase::exp() and most of the time this is the only way it is
398   * used.
399   */
400 template<typename Derived> struct MatrixExponentialReturnValue
401 : public ReturnByValue<MatrixExponentialReturnValue<Derived> >
402 {
403     typedef typename Derived::Index Index;
404   public:
405     /** \brief Constructor.
406       *
407       * \param[in] src %Matrix (expression) forming the argument of the
408       * matrix exponential.
409       */
MatrixExponentialReturnValueMatrixExponentialReturnValue410     MatrixExponentialReturnValue(const Derived& src) : m_src(src) { }
411 
412     /** \brief Compute the matrix exponential.
413       *
414       * \param[out] result the matrix exponential of \p src in the
415       * constructor.
416       */
417     template <typename ResultType>
evalToMatrixExponentialReturnValue418     inline void evalTo(ResultType& result) const
419     {
420       const typename Derived::PlainObject srcEvaluated = m_src.eval();
421       MatrixExponential<typename Derived::PlainObject> me(srcEvaluated);
422       me.compute(result);
423     }
424 
rowsMatrixExponentialReturnValue425     Index rows() const { return m_src.rows(); }
colsMatrixExponentialReturnValue426     Index cols() const { return m_src.cols(); }
427 
428   protected:
429     const Derived& m_src;
430   private:
431     MatrixExponentialReturnValue& operator=(const MatrixExponentialReturnValue&);
432 };
433 
434 namespace internal {
435 template<typename Derived>
436 struct traits<MatrixExponentialReturnValue<Derived> >
437 {
438   typedef typename Derived::PlainObject ReturnType;
439 };
440 }
441 
442 template <typename Derived>
443 const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
444 {
445   eigen_assert(rows() == cols());
446   return MatrixExponentialReturnValue<Derived>(derived());
447 }
448 
449 } // end namespace Eigen
450 
451 #endif // EIGEN_MATRIX_EXPONENTIAL
452