1// This file is part of Eigen, a lightweight C++ template library
4// Copyright (C) 2009 Jitse Niesen <jitse@maths.leeds.ac.uk>
5// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
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/.
11#ifndef EIGEN_MATRIX_FUNCTIONS
12#define EIGEN_MATRIX_FUNCTIONS
17#include "../../Eigen/Core"
18#include "../../Eigen/LU"
19#include "../../Eigen/Eigenvalues"
22 * \defgroup MatrixFunctions_Module Matrix functions module
23 * \brief This module aims to provide various methods for the computation of
26 * To use this module, add
28 * #include <unsupported/Eigen/MatrixFunctions>
30 * at the start of your source file.
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
43 * These methods are the main entry points to this module.
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
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]
56#include "../../Eigen/src/Core/util/DisableStupidWarnings.h"
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"
64#include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"
68\page matrixbaseextra_page
69\ingroup MatrixFunctions_Module
71\section matrixbaseextra MatrixBase methods defined in the MatrixFunctions module
73The remainder of the page documents the following MatrixBase methods
74which are defined in the MatrixFunctions module.
78\subsection matrixbase_cos MatrixBase::cos()
80Compute the matrix cosine.
83const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
86\param[in] M a square matrix.
87\returns expression representing \f$ \cos(M) \f$.
89This function computes the matrix cosine. Use ArrayBase::cos() for computing the entry-wise cosine.
91The implementation calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cos().
93\sa \ref matrixbase_sin "sin()" for an example.
97\subsection matrixbase_cosh MatrixBase::cosh()
99Compute the matrix hyberbolic cosine.
102const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const
105\param[in] M a square matrix.
106\returns expression representing \f$ \cosh(M) \f$
108This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cosh().
110\sa \ref matrixbase_sinh "sinh()" for an example.
114\subsection matrixbase_exp MatrixBase::exp()
116Compute the matrix exponential.
119const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
122\param[in] M matrix whose exponential is to be computed.
123\returns expression representing the matrix exponential of \p M.
125The matrix exponential of \f$ M \f$ is defined by
126\f[ \exp(M) = \sum_{k=0}^\infty \frac{M^k}{k!}. \f]
127The matrix exponential can be used to solve linear ordinary
128differential equations: the solution of \f$ y' = My \f$ with the
129initial condition \f$ y(0) = y_0 \f$ is given by
130\f$ y(t) = \exp(M) y_0 \f$.
132The matrix exponential is different from applying the exp function to all the entries in the matrix.
133Use ArrayBase::exp() if you want to do the latter.
135The cost of the computation is approximately \f$ 20 n^3 \f$ for
136matrices of size \f$ n \f$. The number 20 depends weakly on the
139The matrix exponential is computed using the scaling-and-squaring
140method combined with Padé approximation. The matrix is first
141rescaled, then the exponential of the reduced matrix is computed
142approximant, and then the rescaling is undone by repeated
143squaring. The degree of the Padé approximant is chosen such
144that the approximation error is less than the round-off
145error. However, errors may accumulate during the squaring phase.
147Details of the algorithm can be found in: Nicholas J. Higham, "The
148scaling and squaring method for the matrix exponential revisited,"
149<em>SIAM J. %Matrix Anal. Applic.</em>, <b>26</b>:1179–1193,
152Example: The following program checks that
153\f[ \exp \left[ \begin{array}{ccc}
154 0 & \frac14\pi & 0 \\
155 -\frac14\pi & 0 & 0 \\
157 \end{array} \right] = \left[ \begin{array}{ccc}
158 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
159 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
161 \end{array} \right]. \f]
162This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
165\include MatrixExponential.cpp
166Output: \verbinclude MatrixExponential.out
168\note \p M has to be a matrix of \c float, \c double, `long double`
169\c complex<float>, \c complex<double>, or `complex<long double>` .
172\subsection matrixbase_log MatrixBase::log()
174Compute the matrix logarithm.
177const MatrixLogarithmReturnValue<Derived> MatrixBase<Derived>::log() const
180\param[in] M invertible matrix whose logarithm is to be computed.
181\returns expression representing the matrix logarithm root of \p M.
183The matrix logarithm of \f$ M \f$ is a matrix \f$ X \f$ such that
184\f$ \exp(X) = M \f$ where exp denotes the matrix exponential. As for
185the scalar logarithm, the equation \f$ \exp(X) = M \f$ may have
186multiple solutions; this function returns a matrix whose eigenvalues
187have imaginary part in the interval \f$ (-\pi,\pi] \f$.
189The matrix logarithm is different from applying the log function to all the entries in the matrix.
190Use ArrayBase::log() if you want to do the latter.
192In the real case, the matrix \f$ M \f$ should be invertible and
193it should have no eigenvalues which are real and negative (pairs of
194complex conjugate eigenvalues are allowed). In the complex case, it
195only needs to be invertible.
197This function computes the matrix logarithm using the Schur-Parlett
198algorithm as implemented by MatrixBase::matrixFunction(). The
199logarithm of an atomic block is computed by MatrixLogarithmAtomic,
200which uses direct computation for 1-by-1 and 2-by-2 blocks and an
201inverse scaling-and-squaring algorithm for bigger blocks, with the
202square roots computed by MatrixBase::sqrt().
204Details of the algorithm can be found in Section 11.6.2 of:
206<em>Functions of Matrices: Theory and Computation</em>,
207SIAM 2008. ISBN 978-0-898716-46-7.
209Example: The following program checks that
210\f[ \log \left[ \begin{array}{ccc}
211 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
212 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
214 \end{array} \right] = \left[ \begin{array}{ccc}
215 0 & \frac14\pi & 0 \\
216 -\frac14\pi & 0 & 0 \\
218 \end{array} \right]. \f]
219This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
220the z-axis. This is the inverse of the example used in the
221documentation of \ref matrixbase_exp "exp()".
223\include MatrixLogarithm.cpp
224Output: \verbinclude MatrixLogarithm.out
226\note \p M has to be a matrix of \c float, \c double, `long
227double`, \c complex<float>, \c complex<double>, or `complex<long double>`.
229\sa MatrixBase::exp(), MatrixBase::matrixFunction(),
230 class MatrixLogarithmAtomic, MatrixBase::sqrt().
233\subsection matrixbase_pow MatrixBase::pow()
235Compute the matrix raised to arbitrary real power.
238const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(RealScalar p) const
241\param[in] M base of the matrix power, should be a square matrix.
242\param[in] p exponent of the matrix power.
244The matrix power \f$ M^p \f$ is defined as \f$ \exp(p \log(M)) \f$,
245where exp denotes the matrix exponential, and log denotes the matrix
246logarithm. This is different from raising all the entries in the matrix
247to the p-th power. Use ArrayBase::pow() if you want to do the latter.
249If \p p is complex, the scalar type of \p M should be the type of \p
250p . \f$ M^p \f$ simply evaluates into \f$ \exp(p \log(M)) \f$.
251Therefore, the matrix \f$ M \f$ should meet the conditions to be an
252argument of matrix logarithm.
254If \p p is real, it is casted into the real scalar type of \p M. Then
255this function computes the matrix power using the Schur-Padé
256algorithm as implemented by class MatrixPower. The exponent is split
257into integral part and fractional part, where the fractional part is
258in the interval \f$ (-1, 1) \f$. The main diagonal and the first
259super-diagonal is directly computed.
261If \p M is singular with a semisimple zero eigenvalue and \p p is
262positive, the Schur factor \f$ T \f$ is reordered with Givens
265\f[ T = \left[ \begin{array}{cc}
268 \end{array} \right] \f]
270where \f$ T_1 \f$ is invertible. Then \f$ T^p \f$ is given by
272\f[ T^p = \left[ \begin{array}{cc}
273 T_1^p & T_1^{-1} T_1^p T_2 \\
275 \end{array}. \right] \f]
277\warning Fractional power of a matrix with a non-semisimple zero
278eigenvalue is not well-defined. We introduce an assertion failure
279against inaccurate result, e.g. \code
280#include <unsupported/Eigen/MatrixFunctions>
290 std::cout << A.pow(0.37) << std::endl;
292 // The 1 makes eigenvalue 0 non-semisimple.
293 A.coeffRef(0, 1) = 1;
295 // This fails if EIGEN_NO_DEBUG is undefined.
296 std::cout << A.pow(0.37) << std::endl;
302Details of the algorithm can be found in: Nicholas J. Higham and
303Lijing Lin, "A Schur-Padé algorithm for fractional powers of a
304matrix," <em>SIAM J. %Matrix Anal. Applic.</em>,
305<b>32(3)</b>:1056–1078, 2011.
307Example: The following program checks that
308\f[ \left[ \begin{array}{ccc}
309 \cos1 & -\sin1 & 0 \\
312 \end{array} \right]^{\frac14\pi} = \left[ \begin{array}{ccc}
313 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
314 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
316 \end{array} \right]. \f]
317This corresponds to \f$ \frac14\pi \f$ rotations of 1 radian around
320\include MatrixPower.cpp
321Output: \verbinclude MatrixPower.out
323MatrixBase::pow() is user-friendly. However, there are some
324circumstances under which you should use class MatrixPower directly.
325MatrixPower can save the result of Schur decomposition, so it's
326better for computing various powers for the same matrix.
329\include MatrixPower_optimal.cpp
330Output: \verbinclude MatrixPower_optimal.out
332\note \p M has to be a matrix of \c float, \c double, `long
333double`, \c complex<float>, \c complex<double>, or
334\c complex<long double> .
336\sa MatrixBase::exp(), MatrixBase::log(), class MatrixPower.
339\subsection matrixbase_matrixfunction MatrixBase::matrixFunction()
341Compute a matrix function.
344const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const
347\param[in] M argument of matrix function, should be a square matrix.
348\param[in] f an entire function; \c f(x,n) should compute the n-th
350\returns expression representing \p f applied to \p M.
352Suppose that \p M is a matrix whose entries have type \c Scalar.
353Then, the second argument, \p f, should be a function with prototype
355ComplexScalar f(ComplexScalar, int)
357where \c ComplexScalar = \c std::complex<Scalar> if \c Scalar is
358real (e.g., \c float or \c double) and \c ComplexScalar =
359\c Scalar if \c Scalar is complex. The return value of \c f(x,n)
360should be \f$ f^{(n)}(x) \f$, the n-th derivative of f at x.
362This routine uses the algorithm described in:
363Philip Davies and Nicholas J. Higham,
364"A Schur-Parlett algorithm for computing matrix functions",
365<em>SIAM J. %Matrix Anal. Applic.</em>, <b>25</b>:464–485, 2003.
367The actual work is done by the MatrixFunction class.
369Example: The following program checks that
370\f[ \exp \left[ \begin{array}{ccc}
371 0 & \frac14\pi & 0 \\
372 -\frac14\pi & 0 & 0 \\
374 \end{array} \right] = \left[ \begin{array}{ccc}
375 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
376 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
378 \end{array} \right]. \f]
379This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
380the z-axis. This is the same example as used in the documentation
381of \ref matrixbase_exp "exp()".
383\include MatrixFunction.cpp
384Output: \verbinclude MatrixFunction.out
386Note that the function \c expfn is defined for complex numbers
387\c x, even though the matrix \c A is over the reals. Instead of
388\c expfn, we could also have used StdStemFunctions::exp:
390A.matrixFunction(StdStemFunctions<std::complex<double> >::exp, &B);
395\subsection matrixbase_sin MatrixBase::sin()
397Compute the matrix sine.
400const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
403\param[in] M a square matrix.
404\returns expression representing \f$ \sin(M) \f$.
406This function computes the matrix sine. Use ArrayBase::sin() for computing the entry-wise sine.
408The implementation calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sin().
410Example: \include MatrixSine.cpp
411Output: \verbinclude MatrixSine.out
415\subsection matrixbase_sinh MatrixBase::sinh()
417Compute the matrix hyperbolic sine.
420MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const
423\param[in] M a square matrix.
424\returns expression representing \f$ \sinh(M) \f$
426This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sinh().
428Example: \include MatrixSinh.cpp
429Output: \verbinclude MatrixSinh.out
432\subsection matrixbase_sqrt MatrixBase::sqrt()
434Compute the matrix square root.
437const MatrixSquareRootReturnValue<Derived> MatrixBase<Derived>::sqrt() const
440\param[in] M invertible matrix whose square root is to be computed.
441\returns expression representing the matrix square root of \p M.
443The matrix square root of \f$ M \f$ is the matrix \f$ M^{1/2} \f$
444whose square is the original matrix; so if \f$ S = M^{1/2} \f$ then
445\f$ S^2 = M \f$. This is different from taking the square root of all
446the entries in the matrix; use ArrayBase::sqrt() if you want to do the
449In the <b>real case</b>, the matrix \f$ M \f$ should be invertible and
450it should have no eigenvalues which are real and negative (pairs of
451complex conjugate eigenvalues are allowed). In that case, the matrix
452has a square root which is also real, and this is the square root
453computed by this function.
455The matrix square root is computed by first reducing the matrix to
456quasi-triangular form with the real Schur decomposition. The square
457root of the quasi-triangular matrix can then be computed directly. The
458cost is approximately \f$ 25 n^3 \f$ real flops for the real Schur
459decomposition and \f$ 3\frac13 n^3 \f$ real flops for the remainder
460(though the computation time in practice is likely more than this
463Details of the algorithm can be found in: Nicholas J. Highan,
464"Computing real square roots of a real matrix", <em>Linear Algebra
465Appl.</em>, 88/89:405–430, 1987.
467If the matrix is <b>positive-definite symmetric</b>, then the square
468root is also positive-definite symmetric. In this case, it is best to
469use SelfAdjointEigenSolver::operatorSqrt() to compute it.
471In the <b>complex case</b>, the matrix \f$ M \f$ should be invertible;
472this is a restriction of the algorithm. The square root computed by
473this algorithm is the one whose eigenvalues have an argument in the
474interval \f$ (-\frac12\pi, \frac12\pi] \f$. This is the usual branch
477The computation is the same as in the real case, except that the
478complex Schur decomposition is used to reduce the matrix to a
479triangular matrix. The theoretical cost is the same. Details are in:
480Åke Björck and Sven Hammarling, "A Schur method for the
481square root of a matrix", <em>Linear Algebra Appl.</em>,
48252/53:127–140, 1983.
484Example: The following program checks that the square root of
485\f[ \left[ \begin{array}{cc}
486 \cos(\frac13\pi) & -\sin(\frac13\pi) \\
487 \sin(\frac13\pi) & \cos(\frac13\pi)
488 \end{array} \right], \f]
489corresponding to a rotation over 60 degrees, is a rotation over 30 degrees:
490\f[ \left[ \begin{array}{cc}
491 \cos(\frac16\pi) & -\sin(\frac16\pi) \\
492 \sin(\frac16\pi) & \cos(\frac16\pi)
493 \end{array} \right]. \f]
495\include MatrixSquareRoot.cpp
496Output: \verbinclude MatrixSquareRoot.out
498\sa class RealSchur, class ComplexSchur, class MatrixSquareRoot,
499 SelfAdjointEigenSolver::operatorSqrt().
503#endif // EIGEN_MATRIX_FUNCTIONS