Warning, /include/eigen3/unsupported/Eigen/MatrixFunctions is written in an unsupported language. File is not indexed.
0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2009 Jitse Niesen <jitse@maths.leeds.ac.uk>
0005 // Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_MATRIX_FUNCTIONS
0012 #define EIGEN_MATRIX_FUNCTIONS
0013
0014 #include <cfloat>
0015 #include <list>
0016
0017 #include "../../Eigen/Core"
0018 #include "../../Eigen/LU"
0019 #include "../../Eigen/Eigenvalues"
0020
0021 /**
0022 * \defgroup MatrixFunctions_Module Matrix functions module
0023 * \brief This module aims to provide various methods for the computation of
0024 * matrix functions.
0025 *
0026 * To use this module, add
0027 * \code
0028 * #include <unsupported/Eigen/MatrixFunctions>
0029 * \endcode
0030 * at the start of your source file.
0031 *
0032 * This module defines the following MatrixBase methods.
0033 * - \ref matrixbase_cos "MatrixBase::cos()", for computing the matrix cosine
0034 * - \ref matrixbase_cosh "MatrixBase::cosh()", for computing the matrix hyperbolic cosine
0035 * - \ref matrixbase_exp "MatrixBase::exp()", for computing the matrix exponential
0036 * - \ref matrixbase_log "MatrixBase::log()", for computing the matrix logarithm
0037 * - \ref matrixbase_pow "MatrixBase::pow()", for computing the matrix power
0038 * - \ref matrixbase_matrixfunction "MatrixBase::matrixFunction()", for computing general matrix functions
0039 * - \ref matrixbase_sin "MatrixBase::sin()", for computing the matrix sine
0040 * - \ref matrixbase_sinh "MatrixBase::sinh()", for computing the matrix hyperbolic sine
0041 * - \ref matrixbase_sqrt "MatrixBase::sqrt()", for computing the matrix square root
0042 *
0043 * These methods are the main entry points to this module.
0044 *
0045 * %Matrix functions are defined as follows. Suppose that \f$ f \f$
0046 * is an entire function (that is, a function on the complex plane
0047 * that is everywhere complex differentiable). Then its Taylor
0048 * series
0049 * \f[ f(0) + f'(0) x + \frac{f''(0)}{2} x^2 + \frac{f'''(0)}{3!} x^3 + \cdots \f]
0050 * converges to \f$ f(x) \f$. In this case, we can define the matrix
0051 * function by the same series:
0052 * \f[ f(M) = f(0) + f'(0) M + \frac{f''(0)}{2} M^2 + \frac{f'''(0)}{3!} M^3 + \cdots \f]
0053 *
0054 */
0055
0056 #include "../../Eigen/src/Core/util/DisableStupidWarnings.h"
0057
0058 #include "src/MatrixFunctions/MatrixExponential.h"
0059 #include "src/MatrixFunctions/MatrixFunction.h"
0060 #include "src/MatrixFunctions/MatrixSquareRoot.h"
0061 #include "src/MatrixFunctions/MatrixLogarithm.h"
0062 #include "src/MatrixFunctions/MatrixPower.h"
0063
0064 #include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"
0065
0066
0067 /**
0068 \page matrixbaseextra_page
0069 \ingroup MatrixFunctions_Module
0070
0071 \section matrixbaseextra MatrixBase methods defined in the MatrixFunctions module
0072
0073 The remainder of the page documents the following MatrixBase methods
0074 which are defined in the MatrixFunctions module.
0075
0076
0077
0078 \subsection matrixbase_cos MatrixBase::cos()
0079
0080 Compute the matrix cosine.
0081
0082 \code
0083 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
0084 \endcode
0085
0086 \param[in] M a square matrix.
0087 \returns expression representing \f$ \cos(M) \f$.
0088
0089 This function computes the matrix cosine. Use ArrayBase::cos() for computing the entry-wise cosine.
0090
0091 The implementation calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cos().
0092
0093 \sa \ref matrixbase_sin "sin()" for an example.
0094
0095
0096
0097 \subsection matrixbase_cosh MatrixBase::cosh()
0098
0099 Compute the matrix hyberbolic cosine.
0100
0101 \code
0102 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const
0103 \endcode
0104
0105 \param[in] M a square matrix.
0106 \returns expression representing \f$ \cosh(M) \f$
0107
0108 This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cosh().
0109
0110 \sa \ref matrixbase_sinh "sinh()" for an example.
0111
0112
0113
0114 \subsection matrixbase_exp MatrixBase::exp()
0115
0116 Compute the matrix exponential.
0117
0118 \code
0119 const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
0120 \endcode
0121
0122 \param[in] M matrix whose exponential is to be computed.
0123 \returns expression representing the matrix exponential of \p M.
0124
0125 The matrix exponential of \f$ M \f$ is defined by
0126 \f[ \exp(M) = \sum_{k=0}^\infty \frac{M^k}{k!}. \f]
0127 The matrix exponential can be used to solve linear ordinary
0128 differential equations: the solution of \f$ y' = My \f$ with the
0129 initial condition \f$ y(0) = y_0 \f$ is given by
0130 \f$ y(t) = \exp(M) y_0 \f$.
0131
0132 The matrix exponential is different from applying the exp function to all the entries in the matrix.
0133 Use ArrayBase::exp() if you want to do the latter.
0134
0135 The cost of the computation is approximately \f$ 20 n^3 \f$ for
0136 matrices of size \f$ n \f$. The number 20 depends weakly on the
0137 norm of the matrix.
0138
0139 The matrix exponential is computed using the scaling-and-squaring
0140 method combined with Padé approximation. The matrix is first
0141 rescaled, then the exponential of the reduced matrix is computed
0142 approximant, and then the rescaling is undone by repeated
0143 squaring. The degree of the Padé approximant is chosen such
0144 that the approximation error is less than the round-off
0145 error. However, errors may accumulate during the squaring phase.
0146
0147 Details of the algorithm can be found in: Nicholas J. Higham, "The
0148 scaling and squaring method for the matrix exponential revisited,"
0149 <em>SIAM J. %Matrix Anal. Applic.</em>, <b>26</b>:1179–1193,
0150 2005.
0151
0152 Example: The following program checks that
0153 \f[ \exp \left[ \begin{array}{ccc}
0154 0 & \frac14\pi & 0 \\
0155 -\frac14\pi & 0 & 0 \\
0156 0 & 0 & 0
0157 \end{array} \right] = \left[ \begin{array}{ccc}
0158 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
0159 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
0160 0 & 0 & 1
0161 \end{array} \right]. \f]
0162 This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
0163 the z-axis.
0164
0165 \include MatrixExponential.cpp
0166 Output: \verbinclude MatrixExponential.out
0167
0168 \note \p M has to be a matrix of \c float, \c double, `long double`
0169 \c complex<float>, \c complex<double>, or `complex<long double>` .
0170
0171
0172 \subsection matrixbase_log MatrixBase::log()
0173
0174 Compute the matrix logarithm.
0175
0176 \code
0177 const MatrixLogarithmReturnValue<Derived> MatrixBase<Derived>::log() const
0178 \endcode
0179
0180 \param[in] M invertible matrix whose logarithm is to be computed.
0181 \returns expression representing the matrix logarithm root of \p M.
0182
0183 The matrix logarithm of \f$ M \f$ is a matrix \f$ X \f$ such that
0184 \f$ \exp(X) = M \f$ where exp denotes the matrix exponential. As for
0185 the scalar logarithm, the equation \f$ \exp(X) = M \f$ may have
0186 multiple solutions; this function returns a matrix whose eigenvalues
0187 have imaginary part in the interval \f$ (-\pi,\pi] \f$.
0188
0189 The matrix logarithm is different from applying the log function to all the entries in the matrix.
0190 Use ArrayBase::log() if you want to do the latter.
0191
0192 In the real case, the matrix \f$ M \f$ should be invertible and
0193 it should have no eigenvalues which are real and negative (pairs of
0194 complex conjugate eigenvalues are allowed). In the complex case, it
0195 only needs to be invertible.
0196
0197 This function computes the matrix logarithm using the Schur-Parlett
0198 algorithm as implemented by MatrixBase::matrixFunction(). The
0199 logarithm of an atomic block is computed by MatrixLogarithmAtomic,
0200 which uses direct computation for 1-by-1 and 2-by-2 blocks and an
0201 inverse scaling-and-squaring algorithm for bigger blocks, with the
0202 square roots computed by MatrixBase::sqrt().
0203
0204 Details of the algorithm can be found in Section 11.6.2 of:
0205 Nicholas J. Higham,
0206 <em>Functions of Matrices: Theory and Computation</em>,
0207 SIAM 2008. ISBN 978-0-898716-46-7.
0208
0209 Example: The following program checks that
0210 \f[ \log \left[ \begin{array}{ccc}
0211 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
0212 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
0213 0 & 0 & 1
0214 \end{array} \right] = \left[ \begin{array}{ccc}
0215 0 & \frac14\pi & 0 \\
0216 -\frac14\pi & 0 & 0 \\
0217 0 & 0 & 0
0218 \end{array} \right]. \f]
0219 This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
0220 the z-axis. This is the inverse of the example used in the
0221 documentation of \ref matrixbase_exp "exp()".
0222
0223 \include MatrixLogarithm.cpp
0224 Output: \verbinclude MatrixLogarithm.out
0225
0226 \note \p M has to be a matrix of \c float, \c double, `long
0227 double`, \c complex<float>, \c complex<double>, or `complex<long double>`.
0228
0229 \sa MatrixBase::exp(), MatrixBase::matrixFunction(),
0230 class MatrixLogarithmAtomic, MatrixBase::sqrt().
0231
0232
0233 \subsection matrixbase_pow MatrixBase::pow()
0234
0235 Compute the matrix raised to arbitrary real power.
0236
0237 \code
0238 const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(RealScalar p) const
0239 \endcode
0240
0241 \param[in] M base of the matrix power, should be a square matrix.
0242 \param[in] p exponent of the matrix power.
0243
0244 The matrix power \f$ M^p \f$ is defined as \f$ \exp(p \log(M)) \f$,
0245 where exp denotes the matrix exponential, and log denotes the matrix
0246 logarithm. This is different from raising all the entries in the matrix
0247 to the p-th power. Use ArrayBase::pow() if you want to do the latter.
0248
0249 If \p p is complex, the scalar type of \p M should be the type of \p
0250 p . \f$ M^p \f$ simply evaluates into \f$ \exp(p \log(M)) \f$.
0251 Therefore, the matrix \f$ M \f$ should meet the conditions to be an
0252 argument of matrix logarithm.
0253
0254 If \p p is real, it is casted into the real scalar type of \p M. Then
0255 this function computes the matrix power using the Schur-Padé
0256 algorithm as implemented by class MatrixPower. The exponent is split
0257 into integral part and fractional part, where the fractional part is
0258 in the interval \f$ (-1, 1) \f$. The main diagonal and the first
0259 super-diagonal is directly computed.
0260
0261 If \p M is singular with a semisimple zero eigenvalue and \p p is
0262 positive, the Schur factor \f$ T \f$ is reordered with Givens
0263 rotations, i.e.
0264
0265 \f[ T = \left[ \begin{array}{cc}
0266 T_1 & T_2 \\
0267 0 & 0
0268 \end{array} \right] \f]
0269
0270 where \f$ T_1 \f$ is invertible. Then \f$ T^p \f$ is given by
0271
0272 \f[ T^p = \left[ \begin{array}{cc}
0273 T_1^p & T_1^{-1} T_1^p T_2 \\
0274 0 & 0
0275 \end{array}. \right] \f]
0276
0277 \warning Fractional power of a matrix with a non-semisimple zero
0278 eigenvalue is not well-defined. We introduce an assertion failure
0279 against inaccurate result, e.g. \code
0280 #include <unsupported/Eigen/MatrixFunctions>
0281 #include <iostream>
0282
0283 int main()
0284 {
0285 Eigen::Matrix4d A;
0286 A << 0, 0, 2, 3,
0287 0, 0, 4, 5,
0288 0, 0, 6, 7,
0289 0, 0, 8, 9;
0290 std::cout << A.pow(0.37) << std::endl;
0291
0292 // The 1 makes eigenvalue 0 non-semisimple.
0293 A.coeffRef(0, 1) = 1;
0294
0295 // This fails if EIGEN_NO_DEBUG is undefined.
0296 std::cout << A.pow(0.37) << std::endl;
0297
0298 return 0;
0299 }
0300 \endcode
0301
0302 Details of the algorithm can be found in: Nicholas J. Higham and
0303 Lijing Lin, "A Schur-Padé algorithm for fractional powers of a
0304 matrix," <em>SIAM J. %Matrix Anal. Applic.</em>,
0305 <b>32(3)</b>:1056–1078, 2011.
0306
0307 Example: The following program checks that
0308 \f[ \left[ \begin{array}{ccc}
0309 \cos1 & -\sin1 & 0 \\
0310 \sin1 & \cos1 & 0 \\
0311 0 & 0 & 1
0312 \end{array} \right]^{\frac14\pi} = \left[ \begin{array}{ccc}
0313 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
0314 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
0315 0 & 0 & 1
0316 \end{array} \right]. \f]
0317 This corresponds to \f$ \frac14\pi \f$ rotations of 1 radian around
0318 the z-axis.
0319
0320 \include MatrixPower.cpp
0321 Output: \verbinclude MatrixPower.out
0322
0323 MatrixBase::pow() is user-friendly. However, there are some
0324 circumstances under which you should use class MatrixPower directly.
0325 MatrixPower can save the result of Schur decomposition, so it's
0326 better for computing various powers for the same matrix.
0327
0328 Example:
0329 \include MatrixPower_optimal.cpp
0330 Output: \verbinclude MatrixPower_optimal.out
0331
0332 \note \p M has to be a matrix of \c float, \c double, `long
0333 double`, \c complex<float>, \c complex<double>, or
0334 \c complex<long double> .
0335
0336 \sa MatrixBase::exp(), MatrixBase::log(), class MatrixPower.
0337
0338
0339 \subsection matrixbase_matrixfunction MatrixBase::matrixFunction()
0340
0341 Compute a matrix function.
0342
0343 \code
0344 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const
0345 \endcode
0346
0347 \param[in] M argument of matrix function, should be a square matrix.
0348 \param[in] f an entire function; \c f(x,n) should compute the n-th
0349 derivative of f at x.
0350 \returns expression representing \p f applied to \p M.
0351
0352 Suppose that \p M is a matrix whose entries have type \c Scalar.
0353 Then, the second argument, \p f, should be a function with prototype
0354 \code
0355 ComplexScalar f(ComplexScalar, int)
0356 \endcode
0357 where \c ComplexScalar = \c std::complex<Scalar> if \c Scalar is
0358 real (e.g., \c float or \c double) and \c ComplexScalar =
0359 \c Scalar if \c Scalar is complex. The return value of \c f(x,n)
0360 should be \f$ f^{(n)}(x) \f$, the n-th derivative of f at x.
0361
0362 This routine uses the algorithm described in:
0363 Philip Davies and Nicholas J. Higham,
0364 "A Schur-Parlett algorithm for computing matrix functions",
0365 <em>SIAM J. %Matrix Anal. Applic.</em>, <b>25</b>:464–485, 2003.
0366
0367 The actual work is done by the MatrixFunction class.
0368
0369 Example: The following program checks that
0370 \f[ \exp \left[ \begin{array}{ccc}
0371 0 & \frac14\pi & 0 \\
0372 -\frac14\pi & 0 & 0 \\
0373 0 & 0 & 0
0374 \end{array} \right] = \left[ \begin{array}{ccc}
0375 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
0376 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
0377 0 & 0 & 1
0378 \end{array} \right]. \f]
0379 This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
0380 the z-axis. This is the same example as used in the documentation
0381 of \ref matrixbase_exp "exp()".
0382
0383 \include MatrixFunction.cpp
0384 Output: \verbinclude MatrixFunction.out
0385
0386 Note that the function \c expfn is defined for complex numbers
0387 \c x, even though the matrix \c A is over the reals. Instead of
0388 \c expfn, we could also have used StdStemFunctions::exp:
0389 \code
0390 A.matrixFunction(StdStemFunctions<std::complex<double> >::exp, &B);
0391 \endcode
0392
0393
0394
0395 \subsection matrixbase_sin MatrixBase::sin()
0396
0397 Compute the matrix sine.
0398
0399 \code
0400 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
0401 \endcode
0402
0403 \param[in] M a square matrix.
0404 \returns expression representing \f$ \sin(M) \f$.
0405
0406 This function computes the matrix sine. Use ArrayBase::sin() for computing the entry-wise sine.
0407
0408 The implementation calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sin().
0409
0410 Example: \include MatrixSine.cpp
0411 Output: \verbinclude MatrixSine.out
0412
0413
0414
0415 \subsection matrixbase_sinh MatrixBase::sinh()
0416
0417 Compute the matrix hyperbolic sine.
0418
0419 \code
0420 MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const
0421 \endcode
0422
0423 \param[in] M a square matrix.
0424 \returns expression representing \f$ \sinh(M) \f$
0425
0426 This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sinh().
0427
0428 Example: \include MatrixSinh.cpp
0429 Output: \verbinclude MatrixSinh.out
0430
0431
0432 \subsection matrixbase_sqrt MatrixBase::sqrt()
0433
0434 Compute the matrix square root.
0435
0436 \code
0437 const MatrixSquareRootReturnValue<Derived> MatrixBase<Derived>::sqrt() const
0438 \endcode
0439
0440 \param[in] M invertible matrix whose square root is to be computed.
0441 \returns expression representing the matrix square root of \p M.
0442
0443 The matrix square root of \f$ M \f$ is the matrix \f$ M^{1/2} \f$
0444 whose square is the original matrix; so if \f$ S = M^{1/2} \f$ then
0445 \f$ S^2 = M \f$. This is different from taking the square root of all
0446 the entries in the matrix; use ArrayBase::sqrt() if you want to do the
0447 latter.
0448
0449 In the <b>real case</b>, the matrix \f$ M \f$ should be invertible and
0450 it should have no eigenvalues which are real and negative (pairs of
0451 complex conjugate eigenvalues are allowed). In that case, the matrix
0452 has a square root which is also real, and this is the square root
0453 computed by this function.
0454
0455 The matrix square root is computed by first reducing the matrix to
0456 quasi-triangular form with the real Schur decomposition. The square
0457 root of the quasi-triangular matrix can then be computed directly. The
0458 cost is approximately \f$ 25 n^3 \f$ real flops for the real Schur
0459 decomposition and \f$ 3\frac13 n^3 \f$ real flops for the remainder
0460 (though the computation time in practice is likely more than this
0461 indicates).
0462
0463 Details of the algorithm can be found in: Nicholas J. Highan,
0464 "Computing real square roots of a real matrix", <em>Linear Algebra
0465 Appl.</em>, 88/89:405–430, 1987.
0466
0467 If the matrix is <b>positive-definite symmetric</b>, then the square
0468 root is also positive-definite symmetric. In this case, it is best to
0469 use SelfAdjointEigenSolver::operatorSqrt() to compute it.
0470
0471 In the <b>complex case</b>, the matrix \f$ M \f$ should be invertible;
0472 this is a restriction of the algorithm. The square root computed by
0473 this algorithm is the one whose eigenvalues have an argument in the
0474 interval \f$ (-\frac12\pi, \frac12\pi] \f$. This is the usual branch
0475 cut.
0476
0477 The computation is the same as in the real case, except that the
0478 complex Schur decomposition is used to reduce the matrix to a
0479 triangular matrix. The theoretical cost is the same. Details are in:
0480 Åke Björck and Sven Hammarling, "A Schur method for the
0481 square root of a matrix", <em>Linear Algebra Appl.</em>,
0482 52/53:127–140, 1983.
0483
0484 Example: The following program checks that the square root of
0485 \f[ \left[ \begin{array}{cc}
0486 \cos(\frac13\pi) & -\sin(\frac13\pi) \\
0487 \sin(\frac13\pi) & \cos(\frac13\pi)
0488 \end{array} \right], \f]
0489 corresponding to a rotation over 60 degrees, is a rotation over 30 degrees:
0490 \f[ \left[ \begin{array}{cc}
0491 \cos(\frac16\pi) & -\sin(\frac16\pi) \\
0492 \sin(\frac16\pi) & \cos(\frac16\pi)
0493 \end{array} \right]. \f]
0494
0495 \include MatrixSquareRoot.cpp
0496 Output: \verbinclude MatrixSquareRoot.out
0497
0498 \sa class RealSchur, class ComplexSchur, class MatrixSquareRoot,
0499 SelfAdjointEigenSolver::operatorSqrt().
0500
0501 */
0502
0503 #endif // EIGEN_MATRIX_FUNCTIONS
0504