1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2009 Jitse Niesen <jitse@maths.leeds.ac.uk>
5// Copyright (C) 2012 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_FUNCTIONS
12#define EIGEN_MATRIX_FUNCTIONS
13
14#include <cfloat>
15#include <list>
16
17#include <Eigen/Core>
18#include <Eigen/LU>
19#include <Eigen/Eigenvalues>
20
21/**
22  * \defgroup MatrixFunctions_Module Matrix functions module
23  * \brief This module aims to provide various methods for the computation of
24  * matrix functions.
25  *
26  * To use this module, add
27  * \code
28  * #include <unsupported/Eigen/MatrixFunctions>
29  * \endcode
30  * at the start of your source file.
31  *
32  * This module defines the following MatrixBase methods.
33  *  - \ref matrixbase_cos "MatrixBase::cos()", for computing the matrix cosine
34  *  - \ref matrixbase_cosh "MatrixBase::cosh()", for computing the matrix hyperbolic cosine
35  *  - \ref matrixbase_exp "MatrixBase::exp()", for computing the matrix exponential
36  *  - \ref matrixbase_log "MatrixBase::log()", for computing the matrix logarithm
37  *  - \ref matrixbase_pow "MatrixBase::pow()", for computing the matrix power
38  *  - \ref matrixbase_matrixfunction "MatrixBase::matrixFunction()", for computing general matrix functions
39  *  - \ref matrixbase_sin "MatrixBase::sin()", for computing the matrix sine
40  *  - \ref matrixbase_sinh "MatrixBase::sinh()", for computing the matrix hyperbolic sine
41  *  - \ref matrixbase_sqrt "MatrixBase::sqrt()", for computing the matrix square root
42  *
43  * These methods are the main entry points to this module.
44  *
45  * %Matrix functions are defined as follows.  Suppose that \f$ f \f$
46  * is an entire function (that is, a function on the complex plane
47  * that is everywhere complex differentiable).  Then its Taylor
48  * series
49  * \f[ f(0) + f'(0) x + \frac{f''(0)}{2} x^2 + \frac{f'''(0)}{3!} x^3 + \cdots \f]
50  * converges to \f$ f(x) \f$. In this case, we can define the matrix
51  * function by the same series:
52  * \f[ f(M) = f(0) + f'(0) M + \frac{f''(0)}{2} M^2 + \frac{f'''(0)}{3!} M^3 + \cdots \f]
53  *
54  */
55
56#include "src/MatrixFunctions/MatrixExponential.h"
57#include "src/MatrixFunctions/MatrixFunction.h"
58#include "src/MatrixFunctions/MatrixSquareRoot.h"
59#include "src/MatrixFunctions/MatrixLogarithm.h"
60#include "src/MatrixFunctions/MatrixPower.h"
61
62
63/**
64\page matrixbaseextra_page
65\ingroup MatrixFunctions_Module
66
67\section matrixbaseextra MatrixBase methods defined in the MatrixFunctions module
68
69The remainder of the page documents the following MatrixBase methods
70which are defined in the MatrixFunctions module.
71
72
73
74\subsection matrixbase_cos MatrixBase::cos()
75
76Compute the matrix cosine.
77
78\code
79const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
80\endcode
81
82\param[in]  M  a square matrix.
83\returns  expression representing \f$ \cos(M) \f$.
84
85This function computes the matrix cosine. Use ArrayBase::cos() for computing the entry-wise cosine.
86
87The implementation calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cos().
88
89\sa \ref matrixbase_sin "sin()" for an example.
90
91
92
93\subsection matrixbase_cosh MatrixBase::cosh()
94
95Compute the matrix hyberbolic cosine.
96
97\code
98const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const
99\endcode
100
101\param[in]  M  a square matrix.
102\returns  expression representing \f$ \cosh(M) \f$
103
104This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cosh().
105
106\sa \ref matrixbase_sinh "sinh()" for an example.
107
108
109
110\subsection matrixbase_exp MatrixBase::exp()
111
112Compute the matrix exponential.
113
114\code
115const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
116\endcode
117
118\param[in]  M  matrix whose exponential is to be computed.
119\returns    expression representing the matrix exponential of \p M.
120
121The matrix exponential of \f$ M \f$ is defined by
122\f[ \exp(M) = \sum_{k=0}^\infty \frac{M^k}{k!}. \f]
123The matrix exponential can be used to solve linear ordinary
124differential equations: the solution of \f$ y' = My \f$ with the
125initial condition \f$ y(0) = y_0 \f$ is given by
126\f$ y(t) = \exp(M) y_0 \f$.
127
128The matrix exponential is different from applying the exp function to all the entries in the matrix.
129Use ArrayBase::exp() if you want to do the latter.
130
131The cost of the computation is approximately \f$ 20 n^3 \f$ for
132matrices of size \f$ n \f$. The number 20 depends weakly on the
133norm of the matrix.
134
135The matrix exponential is computed using the scaling-and-squaring
136method combined with Pad&eacute; approximation. The matrix is first
137rescaled, then the exponential of the reduced matrix is computed
138approximant, and then the rescaling is undone by repeated
139squaring. The degree of the Pad&eacute; approximant is chosen such
140that the approximation error is less than the round-off
141error. However, errors may accumulate during the squaring phase.
142
143Details of the algorithm can be found in: Nicholas J. Higham, "The
144scaling and squaring method for the matrix exponential revisited,"
145<em>SIAM J. %Matrix Anal. Applic.</em>, <b>26</b>:1179&ndash;1193,
1462005.
147
148Example: The following program checks that
149\f[ \exp \left[ \begin{array}{ccc}
150      0 & \frac14\pi & 0 \\
151      -\frac14\pi & 0 & 0 \\
152      0 & 0 & 0
153    \end{array} \right] = \left[ \begin{array}{ccc}
154      \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
155      \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
156      0 & 0 & 1
157    \end{array} \right]. \f]
158This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
159the z-axis.
160
161\include MatrixExponential.cpp
162Output: \verbinclude MatrixExponential.out
163
164\note \p M has to be a matrix of \c float, \c double, \c long double
165\c complex<float>, \c complex<double>, or \c complex<long double> .
166
167
168\subsection matrixbase_log MatrixBase::log()
169
170Compute the matrix logarithm.
171
172\code
173const MatrixLogarithmReturnValue<Derived> MatrixBase<Derived>::log() const
174\endcode
175
176\param[in]  M  invertible matrix whose logarithm is to be computed.
177\returns    expression representing the matrix logarithm root of \p M.
178
179The matrix logarithm of \f$ M \f$ is a matrix \f$ X \f$ such that
180\f$ \exp(X) = M \f$ where exp denotes the matrix exponential. As for
181the scalar logarithm, the equation \f$ \exp(X) = M \f$ may have
182multiple solutions; this function returns a matrix whose eigenvalues
183have imaginary part in the interval \f$ (-\pi,\pi] \f$.
184
185The matrix logarithm is different from applying the log function to all the entries in the matrix.
186Use ArrayBase::log() if you want to do the latter.
187
188In the real case, the matrix \f$ M \f$ should be invertible and
189it should have no eigenvalues which are real and negative (pairs of
190complex conjugate eigenvalues are allowed). In the complex case, it
191only needs to be invertible.
192
193This function computes the matrix logarithm using the Schur-Parlett
194algorithm as implemented by MatrixBase::matrixFunction(). The
195logarithm of an atomic block is computed by MatrixLogarithmAtomic,
196which uses direct computation for 1-by-1 and 2-by-2 blocks and an
197inverse scaling-and-squaring algorithm for bigger blocks, with the
198square roots computed by MatrixBase::sqrt().
199
200Details of the algorithm can be found in Section 11.6.2 of:
201Nicholas J. Higham,
202<em>Functions of Matrices: Theory and Computation</em>,
203SIAM 2008. ISBN 978-0-898716-46-7.
204
205Example: The following program checks that
206\f[ \log \left[ \begin{array}{ccc}
207      \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
208      \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
209      0 & 0 & 1
210    \end{array} \right] = \left[ \begin{array}{ccc}
211      0 & \frac14\pi & 0 \\
212      -\frac14\pi & 0 & 0 \\
213      0 & 0 & 0
214    \end{array} \right]. \f]
215This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
216the z-axis. This is the inverse of the example used in the
217documentation of \ref matrixbase_exp "exp()".
218
219\include MatrixLogarithm.cpp
220Output: \verbinclude MatrixLogarithm.out
221
222\note \p M has to be a matrix of \c float, \c double, <tt>long
223double</tt>, \c complex<float>, \c complex<double>, or \c complex<long
224double> .
225
226\sa MatrixBase::exp(), MatrixBase::matrixFunction(),
227    class MatrixLogarithmAtomic, MatrixBase::sqrt().
228
229
230\subsection matrixbase_pow MatrixBase::pow()
231
232Compute the matrix raised to arbitrary real power.
233
234\code
235const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(RealScalar p) const
236\endcode
237
238\param[in]  M  base of the matrix power, should be a square matrix.
239\param[in]  p  exponent of the matrix power.
240
241The matrix power \f$ M^p \f$ is defined as \f$ \exp(p \log(M)) \f$,
242where exp denotes the matrix exponential, and log denotes the matrix
243logarithm. This is different from raising all the entries in the matrix
244to the p-th power. Use ArrayBase::pow() if you want to do the latter.
245
246If \p p is complex, the scalar type of \p M should be the type of \p
247p . \f$ M^p \f$ simply evaluates into \f$ \exp(p \log(M)) \f$.
248Therefore, the matrix \f$ M \f$ should meet the conditions to be an
249argument of matrix logarithm.
250
251If \p p is real, it is casted into the real scalar type of \p M. Then
252this function computes the matrix power using the Schur-Pad&eacute;
253algorithm as implemented by class MatrixPower. The exponent is split
254into integral part and fractional part, where the fractional part is
255in the interval \f$ (-1, 1) \f$. The main diagonal and the first
256super-diagonal is directly computed.
257
258If \p M is singular with a semisimple zero eigenvalue and \p p is
259positive, the Schur factor \f$ T \f$ is reordered with Givens
260rotations, i.e.
261
262\f[ T = \left[ \begin{array}{cc}
263      T_1 & T_2 \\
264      0   & 0
265    \end{array} \right] \f]
266
267where \f$ T_1 \f$ is invertible. Then \f$ T^p \f$ is given by
268
269\f[ T^p = \left[ \begin{array}{cc}
270      T_1^p & T_1^{-1} T_1^p T_2 \\
271      0     & 0
272    \end{array}. \right] \f]
273
274\warning Fractional power of a matrix with a non-semisimple zero
275eigenvalue is not well-defined. We introduce an assertion failure
276against inaccurate result, e.g. \code
277#include <unsupported/Eigen/MatrixFunctions>
278#include <iostream>
279
280int main()
281{
282  Eigen::Matrix4d A;
283  A << 0, 0, 2, 3,
284       0, 0, 4, 5,
285       0, 0, 6, 7,
286       0, 0, 8, 9;
287  std::cout << A.pow(0.37) << std::endl;
288
289  // The 1 makes eigenvalue 0 non-semisimple.
290  A.coeffRef(0, 1) = 1;
291
292  // This fails if EIGEN_NO_DEBUG is undefined.
293  std::cout << A.pow(0.37) << std::endl;
294
295  return 0;
296}
297\endcode
298
299Details of the algorithm can be found in: Nicholas J. Higham and
300Lijing Lin, "A Schur-Pad&eacute; algorithm for fractional powers of a
301matrix," <em>SIAM J. %Matrix Anal. Applic.</em>,
302<b>32(3)</b>:1056&ndash;1078, 2011.
303
304Example: The following program checks that
305\f[ \left[ \begin{array}{ccc}
306      \cos1 & -\sin1 & 0 \\
307      \sin1 & \cos1 & 0 \\
308      0 & 0 & 1
309    \end{array} \right]^{\frac14\pi} = \left[ \begin{array}{ccc}
310      \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
311      \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
312      0 & 0 & 1
313    \end{array} \right]. \f]
314This corresponds to \f$ \frac14\pi \f$ rotations of 1 radian around
315the z-axis.
316
317\include MatrixPower.cpp
318Output: \verbinclude MatrixPower.out
319
320MatrixBase::pow() is user-friendly. However, there are some
321circumstances under which you should use class MatrixPower directly.
322MatrixPower can save the result of Schur decomposition, so it's
323better for computing various powers for the same matrix.
324
325Example:
326\include MatrixPower_optimal.cpp
327Output: \verbinclude MatrixPower_optimal.out
328
329\note \p M has to be a matrix of \c float, \c double, <tt>long
330double</tt>, \c complex<float>, \c complex<double>, or \c complex<long
331double> .
332
333\sa MatrixBase::exp(), MatrixBase::log(), class MatrixPower.
334
335
336\subsection matrixbase_matrixfunction MatrixBase::matrixFunction()
337
338Compute a matrix function.
339
340\code
341const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const
342\endcode
343
344\param[in]  M  argument of matrix function, should be a square matrix.
345\param[in]  f  an entire function; \c f(x,n) should compute the n-th
346derivative of f at x.
347\returns  expression representing \p f applied to \p M.
348
349Suppose that \p M is a matrix whose entries have type \c Scalar.
350Then, the second argument, \p f, should be a function with prototype
351\code
352ComplexScalar f(ComplexScalar, int)
353\endcode
354where \c ComplexScalar = \c std::complex<Scalar> if \c Scalar is
355real (e.g., \c float or \c double) and \c ComplexScalar =
356\c Scalar if \c Scalar is complex. The return value of \c f(x,n)
357should be \f$ f^{(n)}(x) \f$, the n-th derivative of f at x.
358
359This routine uses the algorithm described in:
360Philip Davies and Nicholas J. Higham,
361"A Schur-Parlett algorithm for computing matrix functions",
362<em>SIAM J. %Matrix Anal. Applic.</em>, <b>25</b>:464&ndash;485, 2003.
363
364The actual work is done by the MatrixFunction class.
365
366Example: The following program checks that
367\f[ \exp \left[ \begin{array}{ccc}
368      0 & \frac14\pi & 0 \\
369      -\frac14\pi & 0 & 0 \\
370      0 & 0 & 0
371    \end{array} \right] = \left[ \begin{array}{ccc}
372      \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
373      \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
374      0 & 0 & 1
375    \end{array} \right]. \f]
376This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
377the z-axis. This is the same example as used in the documentation
378of \ref matrixbase_exp "exp()".
379
380\include MatrixFunction.cpp
381Output: \verbinclude MatrixFunction.out
382
383Note that the function \c expfn is defined for complex numbers
384\c x, even though the matrix \c A is over the reals. Instead of
385\c expfn, we could also have used StdStemFunctions::exp:
386\code
387A.matrixFunction(StdStemFunctions<std::complex<double> >::exp, &B);
388\endcode
389
390
391
392\subsection matrixbase_sin MatrixBase::sin()
393
394Compute the matrix sine.
395
396\code
397const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
398\endcode
399
400\param[in]  M  a square matrix.
401\returns  expression representing \f$ \sin(M) \f$.
402
403This function computes the matrix sine. Use ArrayBase::sin() for computing the entry-wise sine.
404
405The implementation calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sin().
406
407Example: \include MatrixSine.cpp
408Output: \verbinclude MatrixSine.out
409
410
411
412\subsection matrixbase_sinh MatrixBase::sinh()
413
414Compute the matrix hyperbolic sine.
415
416\code
417MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const
418\endcode
419
420\param[in]  M  a square matrix.
421\returns  expression representing \f$ \sinh(M) \f$
422
423This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sinh().
424
425Example: \include MatrixSinh.cpp
426Output: \verbinclude MatrixSinh.out
427
428
429\subsection matrixbase_sqrt MatrixBase::sqrt()
430
431Compute the matrix square root.
432
433\code
434const MatrixSquareRootReturnValue<Derived> MatrixBase<Derived>::sqrt() const
435\endcode
436
437\param[in]  M  invertible matrix whose square root is to be computed.
438\returns    expression representing the matrix square root of \p M.
439
440The matrix square root of \f$ M \f$ is the matrix \f$ M^{1/2} \f$
441whose square is the original matrix; so if \f$ S = M^{1/2} \f$ then
442\f$ S^2 = M \f$. This is different from taking the square root of all
443the entries in the matrix; use ArrayBase::sqrt() if you want to do the
444latter.
445
446In the <b>real case</b>, the matrix \f$ M \f$ should be invertible and
447it should have no eigenvalues which are real and negative (pairs of
448complex conjugate eigenvalues are allowed). In that case, the matrix
449has a square root which is also real, and this is the square root
450computed by this function.
451
452The matrix square root is computed by first reducing the matrix to
453quasi-triangular form with the real Schur decomposition. The square
454root of the quasi-triangular matrix can then be computed directly. The
455cost is approximately \f$ 25 n^3 \f$ real flops for the real Schur
456decomposition and \f$ 3\frac13 n^3 \f$ real flops for the remainder
457(though the computation time in practice is likely more than this
458indicates).
459
460Details of the algorithm can be found in: Nicholas J. Highan,
461"Computing real square roots of a real matrix", <em>Linear Algebra
462Appl.</em>, 88/89:405&ndash;430, 1987.
463
464If the matrix is <b>positive-definite symmetric</b>, then the square
465root is also positive-definite symmetric. In this case, it is best to
466use SelfAdjointEigenSolver::operatorSqrt() to compute it.
467
468In the <b>complex case</b>, the matrix \f$ M \f$ should be invertible;
469this is a restriction of the algorithm. The square root computed by
470this algorithm is the one whose eigenvalues have an argument in the
471interval \f$ (-\frac12\pi, \frac12\pi] \f$. This is the usual branch
472cut.
473
474The computation is the same as in the real case, except that the
475complex Schur decomposition is used to reduce the matrix to a
476triangular matrix. The theoretical cost is the same. Details are in:
477&Aring;ke Bj&ouml;rck and Sven Hammarling, "A Schur method for the
478square root of a matrix", <em>Linear Algebra Appl.</em>,
47952/53:127&ndash;140, 1983.
480
481Example: The following program checks that the square root of
482\f[ \left[ \begin{array}{cc}
483              \cos(\frac13\pi) & -\sin(\frac13\pi) \\
484              \sin(\frac13\pi) & \cos(\frac13\pi)
485    \end{array} \right], \f]
486corresponding to a rotation over 60 degrees, is a rotation over 30 degrees:
487\f[ \left[ \begin{array}{cc}
488              \cos(\frac16\pi) & -\sin(\frac16\pi) \\
489              \sin(\frac16\pi) & \cos(\frac16\pi)
490    \end{array} \right]. \f]
491
492\include MatrixSquareRoot.cpp
493Output: \verbinclude MatrixSquareRoot.out
494
495\sa class RealSchur, class ComplexSchur, class MatrixSquareRoot,
496    SelfAdjointEigenSolver::operatorSqrt().
497
498*/
499
500#endif // EIGEN_MATRIX_FUNCTIONS
501
502