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