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// 6// This Source Code Form is subject to the terms of the Mozilla 7// Public License v. 2.0. If a copy of the MPL was not distributed 8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10#ifndef EIGEN_MATRIX_FUNCTIONS 11#define EIGEN_MATRIX_FUNCTIONS 12 13#include <cfloat> 14#include <list> 15#include <functional> 16#include <iterator> 17 18#include <Eigen/Core> 19#include <Eigen/LU> 20#include <Eigen/Eigenvalues> 21 22/** \ingroup Unsupported_modules 23 * \defgroup MatrixFunctions_Module Matrix functions module 24 * \brief This module aims to provide various methods for the computation of 25 * matrix functions. 26 * 27 * To use this module, add 28 * \code 29 * #include <unsupported/Eigen/MatrixFunctions> 30 * \endcode 31 * at the start of your source file. 32 * 33 * This module defines the following MatrixBase methods. 34 * - \ref matrixbase_cos "MatrixBase::cos()", for computing the matrix cosine 35 * - \ref matrixbase_cosh "MatrixBase::cosh()", for computing the matrix hyperbolic cosine 36 * - \ref matrixbase_exp "MatrixBase::exp()", for computing the matrix exponential 37 * - \ref matrixbase_log "MatrixBase::log()", for computing the matrix logarithm 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 61 62 63/** 64\page matrixbaseextra MatrixBase methods defined in the MatrixFunctions module 65\ingroup MatrixFunctions_Module 66 67The remainder of the page documents the following MatrixBase methods 68which are defined in the MatrixFunctions module. 69 70 71 72\section matrixbase_cos MatrixBase::cos() 73 74Compute the matrix cosine. 75 76\code 77const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const 78\endcode 79 80\param[in] M a square matrix. 81\returns expression representing \f$ \cos(M) \f$. 82 83This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cos(). 84 85\sa \ref matrixbase_sin "sin()" for an example. 86 87 88 89\section matrixbase_cosh MatrixBase::cosh() 90 91Compute the matrix hyberbolic cosine. 92 93\code 94const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const 95\endcode 96 97\param[in] M a square matrix. 98\returns expression representing \f$ \cosh(M) \f$ 99 100This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cosh(). 101 102\sa \ref matrixbase_sinh "sinh()" for an example. 103 104 105 106\section matrixbase_exp MatrixBase::exp() 107 108Compute the matrix exponential. 109 110\code 111const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const 112\endcode 113 114\param[in] M matrix whose exponential is to be computed. 115\returns expression representing the matrix exponential of \p M. 116 117The matrix exponential of \f$ M \f$ is defined by 118\f[ \exp(M) = \sum_{k=0}^\infty \frac{M^k}{k!}. \f] 119The matrix exponential can be used to solve linear ordinary 120differential equations: the solution of \f$ y' = My \f$ with the 121initial condition \f$ y(0) = y_0 \f$ is given by 122\f$ y(t) = \exp(M) y_0 \f$. 123 124The cost of the computation is approximately \f$ 20 n^3 \f$ for 125matrices of size \f$ n \f$. The number 20 depends weakly on the 126norm of the matrix. 127 128The matrix exponential is computed using the scaling-and-squaring 129method combined with Padé approximation. The matrix is first 130rescaled, then the exponential of the reduced matrix is computed 131approximant, and then the rescaling is undone by repeated 132squaring. The degree of the Padé approximant is chosen such 133that the approximation error is less than the round-off 134error. However, errors may accumulate during the squaring phase. 135 136Details of the algorithm can be found in: Nicholas J. Higham, "The 137scaling and squaring method for the matrix exponential revisited," 138<em>SIAM J. %Matrix Anal. Applic.</em>, <b>26</b>:1179–1193, 1392005. 140 141Example: The following program checks that 142\f[ \exp \left[ \begin{array}{ccc} 143 0 & \frac14\pi & 0 \\ 144 -\frac14\pi & 0 & 0 \\ 145 0 & 0 & 0 146 \end{array} \right] = \left[ \begin{array}{ccc} 147 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ 148 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ 149 0 & 0 & 1 150 \end{array} \right]. \f] 151This corresponds to a rotation of \f$ \frac14\pi \f$ radians around 152the z-axis. 153 154\include MatrixExponential.cpp 155Output: \verbinclude MatrixExponential.out 156 157\note \p M has to be a matrix of \c float, \c double, \c long double 158\c complex<float>, \c complex<double>, or \c complex<long double> . 159 160 161\section matrixbase_log MatrixBase::log() 162 163Compute the matrix logarithm. 164 165\code 166const MatrixLogarithmReturnValue<Derived> MatrixBase<Derived>::log() const 167\endcode 168 169\param[in] M invertible matrix whose logarithm is to be computed. 170\returns expression representing the matrix logarithm root of \p M. 171 172The matrix logarithm of \f$ M \f$ is a matrix \f$ X \f$ such that 173\f$ \exp(X) = M \f$ where exp denotes the matrix exponential. As for 174the scalar logarithm, the equation \f$ \exp(X) = M \f$ may have 175multiple solutions; this function returns a matrix whose eigenvalues 176have imaginary part in the interval \f$ (-\pi,\pi] \f$. 177 178In the real case, the matrix \f$ M \f$ should be invertible and 179it should have no eigenvalues which are real and negative (pairs of 180complex conjugate eigenvalues are allowed). In the complex case, it 181only needs to be invertible. 182 183This function computes the matrix logarithm using the Schur-Parlett 184algorithm as implemented by MatrixBase::matrixFunction(). The 185logarithm of an atomic block is computed by MatrixLogarithmAtomic, 186which uses direct computation for 1-by-1 and 2-by-2 blocks and an 187inverse scaling-and-squaring algorithm for bigger blocks, with the 188square roots computed by MatrixBase::sqrt(). 189 190Details of the algorithm can be found in Section 11.6.2 of: 191Nicholas J. Higham, 192<em>Functions of Matrices: Theory and Computation</em>, 193SIAM 2008. ISBN 978-0-898716-46-7. 194 195Example: The following program checks that 196\f[ \log \left[ \begin{array}{ccc} 197 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ 198 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ 199 0 & 0 & 1 200 \end{array} \right] = \left[ \begin{array}{ccc} 201 0 & \frac14\pi & 0 \\ 202 -\frac14\pi & 0 & 0 \\ 203 0 & 0 & 0 204 \end{array} \right]. \f] 205This corresponds to a rotation of \f$ \frac14\pi \f$ radians around 206the z-axis. This is the inverse of the example used in the 207documentation of \ref matrixbase_exp "exp()". 208 209\include MatrixLogarithm.cpp 210Output: \verbinclude MatrixLogarithm.out 211 212\note \p M has to be a matrix of \c float, \c double, \c long double 213\c complex<float>, \c complex<double>, or \c complex<long double> . 214 215\sa MatrixBase::exp(), MatrixBase::matrixFunction(), 216 class MatrixLogarithmAtomic, MatrixBase::sqrt(). 217 218 219\section matrixbase_matrixfunction MatrixBase::matrixFunction() 220 221Compute a matrix function. 222 223\code 224const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const 225\endcode 226 227\param[in] M argument of matrix function, should be a square matrix. 228\param[in] f an entire function; \c f(x,n) should compute the n-th 229derivative of f at x. 230\returns expression representing \p f applied to \p M. 231 232Suppose that \p M is a matrix whose entries have type \c Scalar. 233Then, the second argument, \p f, should be a function with prototype 234\code 235ComplexScalar f(ComplexScalar, int) 236\endcode 237where \c ComplexScalar = \c std::complex<Scalar> if \c Scalar is 238real (e.g., \c float or \c double) and \c ComplexScalar = 239\c Scalar if \c Scalar is complex. The return value of \c f(x,n) 240should be \f$ f^{(n)}(x) \f$, the n-th derivative of f at x. 241 242This routine uses the algorithm described in: 243Philip Davies and Nicholas J. Higham, 244"A Schur-Parlett algorithm for computing matrix functions", 245<em>SIAM J. %Matrix Anal. Applic.</em>, <b>25</b>:464–485, 2003. 246 247The actual work is done by the MatrixFunction class. 248 249Example: The following program checks that 250\f[ \exp \left[ \begin{array}{ccc} 251 0 & \frac14\pi & 0 \\ 252 -\frac14\pi & 0 & 0 \\ 253 0 & 0 & 0 254 \end{array} \right] = \left[ \begin{array}{ccc} 255 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ 256 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ 257 0 & 0 & 1 258 \end{array} \right]. \f] 259This corresponds to a rotation of \f$ \frac14\pi \f$ radians around 260the z-axis. This is the same example as used in the documentation 261of \ref matrixbase_exp "exp()". 262 263\include MatrixFunction.cpp 264Output: \verbinclude MatrixFunction.out 265 266Note that the function \c expfn is defined for complex numbers 267\c x, even though the matrix \c A is over the reals. Instead of 268\c expfn, we could also have used StdStemFunctions::exp: 269\code 270A.matrixFunction(StdStemFunctions<std::complex<double> >::exp, &B); 271\endcode 272 273 274 275\section matrixbase_sin MatrixBase::sin() 276 277Compute the matrix sine. 278 279\code 280const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const 281\endcode 282 283\param[in] M a square matrix. 284\returns expression representing \f$ \sin(M) \f$. 285 286This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sin(). 287 288Example: \include MatrixSine.cpp 289Output: \verbinclude MatrixSine.out 290 291 292 293\section matrixbase_sinh MatrixBase::sinh() 294 295Compute the matrix hyperbolic sine. 296 297\code 298MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const 299\endcode 300 301\param[in] M a square matrix. 302\returns expression representing \f$ \sinh(M) \f$ 303 304This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sinh(). 305 306Example: \include MatrixSinh.cpp 307Output: \verbinclude MatrixSinh.out 308 309 310\section matrixbase_sqrt MatrixBase::sqrt() 311 312Compute the matrix square root. 313 314\code 315const MatrixSquareRootReturnValue<Derived> MatrixBase<Derived>::sqrt() const 316\endcode 317 318\param[in] M invertible matrix whose square root is to be computed. 319\returns expression representing the matrix square root of \p M. 320 321The matrix square root of \f$ M \f$ is the matrix \f$ M^{1/2} \f$ 322whose square is the original matrix; so if \f$ S = M^{1/2} \f$ then 323\f$ S^2 = M \f$. 324 325In the <b>real case</b>, the matrix \f$ M \f$ should be invertible and 326it should have no eigenvalues which are real and negative (pairs of 327complex conjugate eigenvalues are allowed). In that case, the matrix 328has a square root which is also real, and this is the square root 329computed by this function. 330 331The matrix square root is computed by first reducing the matrix to 332quasi-triangular form with the real Schur decomposition. The square 333root of the quasi-triangular matrix can then be computed directly. The 334cost is approximately \f$ 25 n^3 \f$ real flops for the real Schur 335decomposition and \f$ 3\frac13 n^3 \f$ real flops for the remainder 336(though the computation time in practice is likely more than this 337indicates). 338 339Details of the algorithm can be found in: Nicholas J. Highan, 340"Computing real square roots of a real matrix", <em>Linear Algebra 341Appl.</em>, 88/89:405–430, 1987. 342 343If the matrix is <b>positive-definite symmetric</b>, then the square 344root is also positive-definite symmetric. In this case, it is best to 345use SelfAdjointEigenSolver::operatorSqrt() to compute it. 346 347In the <b>complex case</b>, the matrix \f$ M \f$ should be invertible; 348this is a restriction of the algorithm. The square root computed by 349this algorithm is the one whose eigenvalues have an argument in the 350interval \f$ (-\frac12\pi, \frac12\pi] \f$. This is the usual branch 351cut. 352 353The computation is the same as in the real case, except that the 354complex Schur decomposition is used to reduce the matrix to a 355triangular matrix. The theoretical cost is the same. Details are in: 356Åke Björck and Sven Hammarling, "A Schur method for the 357square root of a matrix", <em>Linear Algebra Appl.</em>, 35852/53:127–140, 1983. 359 360Example: The following program checks that the square root of 361\f[ \left[ \begin{array}{cc} 362 \cos(\frac13\pi) & -\sin(\frac13\pi) \\ 363 \sin(\frac13\pi) & \cos(\frac13\pi) 364 \end{array} \right], \f] 365corresponding to a rotation over 60 degrees, is a rotation over 30 degrees: 366\f[ \left[ \begin{array}{cc} 367 \cos(\frac16\pi) & -\sin(\frac16\pi) \\ 368 \sin(\frac16\pi) & \cos(\frac16\pi) 369 \end{array} \right]. \f] 370 371\include MatrixSquareRoot.cpp 372Output: \verbinclude MatrixSquareRoot.out 373 374\sa class RealSchur, class ComplexSchur, class MatrixSquareRoot, 375 SelfAdjointEigenSolver::operatorSqrt(). 376 377*/ 378 379#endif // EIGEN_MATRIX_FUNCTIONS 380 381