Back to home page

EIC code displayed by LXR

 
 

    


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&eacute; 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&eacute; 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&ndash;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&eacute;
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&eacute; algorithm for fractional powers of a
0304 matrix," <em>SIAM J. %Matrix Anal. Applic.</em>,
0305 <b>32(3)</b>:1056&ndash;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&ndash;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&ndash;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 &Aring;ke Bj&ouml;rck and Sven Hammarling, "A Schur method for the
0481 square root of a matrix", <em>Linear Algebra Appl.</em>,
0482 52/53:127&ndash;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