1// This file is part of Eigen, a lightweight C++ template library 2// for linear algebra. 3// 4// Copyright (C) 2009 Jitse Niesen <[email protected]> 5// Copyright (C) 2012 Chen-Pang He <[email protected]> 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 17#include "../../Eigen/Core" 18#include "../../Eigen/LU" 19#include "../../Eigen/Eigenvalues" 20 21/** 22 * \defgroup MatrixFunctions_Module Matrix functions module 23 * \brief This module aims to provide various methods for the computation of 24 * matrix functions. 25 * 26 * To use this module, add 27 * \code 28 * #include <unsupported/Eigen/MatrixFunctions> 29 * \endcode 30 * at the start of your source file. 31 * 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 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 "../../Eigen/src/Core/util/DisableStupidWarnings.h" 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#include "../../Eigen/src/Core/util/ReenableStupidWarnings.h" 65 66 67/** 68\page matrixbaseextra_page 69\ingroup MatrixFunctions_Module 70 71\section matrixbaseextra MatrixBase methods defined in the MatrixFunctions module 72 73The remainder of the page documents the following MatrixBase methods 74which are defined in the MatrixFunctions module. 75 76 77 78\subsection matrixbase_cos MatrixBase::cos() 79 80Compute the matrix cosine. 81 82\code 83const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const 84\endcode 85 86\param[in] M a square matrix. 87\returns expression representing \f$ \cos(M) \f$. 88 89This function computes the matrix cosine. Use ArrayBase::cos() for computing the entry-wise cosine. 90 91The implementation calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cos(). 92 93\sa \ref matrixbase_sin "sin()" for an example. 94 95 96 97\subsection matrixbase_cosh MatrixBase::cosh() 98 99Compute the matrix hyberbolic cosine. 100 101\code 102const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const 103\endcode 104 105\param[in] M a square matrix. 106\returns expression representing \f$ \cosh(M) \f$ 107 108This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cosh(). 109 110\sa \ref matrixbase_sinh "sinh()" for an example. 111 112 113 114\subsection matrixbase_exp MatrixBase::exp() 115 116Compute the matrix exponential. 117 118\code 119const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const 120\endcode 121 122\param[in] M matrix whose exponential is to be computed. 123\returns expression representing the matrix exponential of \p M. 124 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$. 131 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. 134 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 137norm of the matrix. 138 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. 146 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, 1502005. 151 152Example: The following program checks that 153\f[ \exp \left[ \begin{array}{ccc} 154 0 & \frac14\pi & 0 \\ 155 -\frac14\pi & 0 & 0 \\ 156 0 & 0 & 0 157 \end{array} \right] = \left[ \begin{array}{ccc} 158 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ 159 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ 160 0 & 0 & 1 161 \end{array} \right]. \f] 162This corresponds to a rotation of \f$ \frac14\pi \f$ radians around 163the z-axis. 164 165\include MatrixExponential.cpp 166Output: \verbinclude MatrixExponential.out 167 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>` . 170 171 172\subsection matrixbase_log MatrixBase::log() 173 174Compute the matrix logarithm. 175 176\code 177const MatrixLogarithmReturnValue<Derived> MatrixBase<Derived>::log() const 178\endcode 179 180\param[in] M invertible matrix whose logarithm is to be computed. 181\returns expression representing the matrix logarithm root of \p M. 182 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$. 188 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. 191 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. 196 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(). 203 204Details of the algorithm can be found in Section 11.6.2 of: 205Nicholas J. Higham, 206<em>Functions of Matrices: Theory and Computation</em>, 207SIAM 2008. ISBN 978-0-898716-46-7. 208 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 \\ 213 0 & 0 & 1 214 \end{array} \right] = \left[ \begin{array}{ccc} 215 0 & \frac14\pi & 0 \\ 216 -\frac14\pi & 0 & 0 \\ 217 0 & 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()". 222 223\include MatrixLogarithm.cpp 224Output: \verbinclude MatrixLogarithm.out 225 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>`. 228 229\sa MatrixBase::exp(), MatrixBase::matrixFunction(), 230 class MatrixLogarithmAtomic, MatrixBase::sqrt(). 231 232 233\subsection matrixbase_pow MatrixBase::pow() 234 235Compute the matrix raised to arbitrary real power. 236 237\code 238const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(RealScalar p) const 239\endcode 240 241\param[in] M base of the matrix power, should be a square matrix. 242\param[in] p exponent of the matrix power. 243 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. 248 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. 253 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. 260 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 263rotations, i.e. 264 265\f[ T = \left[ \begin{array}{cc} 266 T_1 & T_2 \\ 267 0 & 0 268 \end{array} \right] \f] 269 270where \f$ T_1 \f$ is invertible. Then \f$ T^p \f$ is given by 271 272\f[ T^p = \left[ \begin{array}{cc} 273 T_1^p & T_1^{-1} T_1^p T_2 \\ 274 0 & 0 275 \end{array}. \right] \f] 276 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> 281#include <iostream> 282 283int main() 284{ 285 Eigen::Matrix4d A; 286 A << 0, 0, 2, 3, 287 0, 0, 4, 5, 288 0, 0, 6, 7, 289 0, 0, 8, 9; 290 std::cout << A.pow(0.37) << std::endl; 291 292 // The 1 makes eigenvalue 0 non-semisimple. 293 A.coeffRef(0, 1) = 1; 294 295 // This fails if EIGEN_NO_DEBUG is undefined. 296 std::cout << A.pow(0.37) << std::endl; 297 298 return 0; 299} 300\endcode 301 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. 306 307Example: The following program checks that 308\f[ \left[ \begin{array}{ccc} 309 \cos1 & -\sin1 & 0 \\ 310 \sin1 & \cos1 & 0 \\ 311 0 & 0 & 1 312 \end{array} \right]^{\frac14\pi} = \left[ \begin{array}{ccc} 313 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ 314 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ 315 0 & 0 & 1 316 \end{array} \right]. \f] 317This corresponds to \f$ \frac14\pi \f$ rotations of 1 radian around 318the z-axis. 319 320\include MatrixPower.cpp 321Output: \verbinclude MatrixPower.out 322 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. 327 328Example: 329\include MatrixPower_optimal.cpp 330Output: \verbinclude MatrixPower_optimal.out 331 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> . 335 336\sa MatrixBase::exp(), MatrixBase::log(), class MatrixPower. 337 338 339\subsection matrixbase_matrixfunction MatrixBase::matrixFunction() 340 341Compute a matrix function. 342 343\code 344const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const 345\endcode 346 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 349derivative of f at x. 350\returns expression representing \p f applied to \p M. 351 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 354\code 355ComplexScalar f(ComplexScalar, int) 356\endcode 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. 361 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. 366 367The actual work is done by the MatrixFunction class. 368 369Example: The following program checks that 370\f[ \exp \left[ \begin{array}{ccc} 371 0 & \frac14\pi & 0 \\ 372 -\frac14\pi & 0 & 0 \\ 373 0 & 0 & 0 374 \end{array} \right] = \left[ \begin{array}{ccc} 375 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ 376 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ 377 0 & 0 & 1 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()". 382 383\include MatrixFunction.cpp 384Output: \verbinclude MatrixFunction.out 385 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: 389\code 390A.matrixFunction(StdStemFunctions<std::complex<double> >::exp, &B); 391\endcode 392 393 394 395\subsection matrixbase_sin MatrixBase::sin() 396 397Compute the matrix sine. 398 399\code 400const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const 401\endcode 402 403\param[in] M a square matrix. 404\returns expression representing \f$ \sin(M) \f$. 405 406This function computes the matrix sine. Use ArrayBase::sin() for computing the entry-wise sine. 407 408The implementation calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sin(). 409 410Example: \include MatrixSine.cpp 411Output: \verbinclude MatrixSine.out 412 413 414 415\subsection matrixbase_sinh MatrixBase::sinh() 416 417Compute the matrix hyperbolic sine. 418 419\code 420MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const 421\endcode 422 423\param[in] M a square matrix. 424\returns expression representing \f$ \sinh(M) \f$ 425 426This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sinh(). 427 428Example: \include MatrixSinh.cpp 429Output: \verbinclude MatrixSinh.out 430 431 432\subsection matrixbase_sqrt MatrixBase::sqrt() 433 434Compute the matrix square root. 435 436\code 437const MatrixSquareRootReturnValue<Derived> MatrixBase<Derived>::sqrt() const 438\endcode 439 440\param[in] M invertible matrix whose square root is to be computed. 441\returns expression representing the matrix square root of \p M. 442 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 447latter. 448 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. 454 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 461indicates). 462 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. 466 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. 470 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 475cut. 476 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. 483 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] 494 495\include MatrixSquareRoot.cpp 496Output: \verbinclude MatrixSquareRoot.out 497 498\sa class RealSchur, class ComplexSchur, class MatrixSquareRoot, 499 SelfAdjointEigenSolver::operatorSqrt(). 500 501*/ 502 503#endif // EIGEN_MATRIX_FUNCTIONS 504 505