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// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net> 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#include <functional> 17#include <iterator> 18 19#include <Eigen/Core> 20#include <Eigen/LU> 21#include <Eigen/Eigenvalues> 22 23/** 24 * \defgroup MatrixFunctions_Module Matrix functions module 25 * \brief This module aims to provide various methods for the computation of 26 * matrix functions. 27 * 28 * To use this module, add 29 * \code 30 * #include <unsupported/Eigen/MatrixFunctions> 31 * \endcode 32 * at the start of your source file. 33 * 34 * This module defines the following MatrixBase methods. 35 * - \ref matrixbase_cos "MatrixBase::cos()", for computing the matrix cosine 36 * - \ref matrixbase_cosh "MatrixBase::cosh()", for computing the matrix hyperbolic cosine 37 * - \ref matrixbase_exp "MatrixBase::exp()", for computing the matrix exponential 38 * - \ref matrixbase_log "MatrixBase::log()", for computing the matrix logarithm 39 * - \ref matrixbase_pow "MatrixBase::pow()", for computing the matrix power 40 * - \ref matrixbase_matrixfunction "MatrixBase::matrixFunction()", for computing general matrix functions 41 * - \ref matrixbase_sin "MatrixBase::sin()", for computing the matrix sine 42 * - \ref matrixbase_sinh "MatrixBase::sinh()", for computing the matrix hyperbolic sine 43 * - \ref matrixbase_sqrt "MatrixBase::sqrt()", for computing the matrix square root 44 * 45 * These methods are the main entry points to this module. 46 * 47 * %Matrix functions are defined as follows. Suppose that \f$ f \f$ 48 * is an entire function (that is, a function on the complex plane 49 * that is everywhere complex differentiable). Then its Taylor 50 * series 51 * \f[ f(0) + f'(0) x + \frac{f''(0)}{2} x^2 + \frac{f'''(0)}{3!} x^3 + \cdots \f] 52 * converges to \f$ f(x) \f$. In this case, we can define the matrix 53 * function by the same series: 54 * \f[ f(M) = f(0) + f'(0) M + \frac{f''(0)}{2} M^2 + \frac{f'''(0)}{3!} M^3 + \cdots \f] 55 * 56 */ 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 65/** 66\page matrixbaseextra_page 67\ingroup MatrixFunctions_Module 68 69\section matrixbaseextra MatrixBase methods defined in the MatrixFunctions module 70 71The remainder of the page documents the following MatrixBase methods 72which are defined in the MatrixFunctions module. 73 74 75 76\subsection matrixbase_cos MatrixBase::cos() 77 78Compute the matrix cosine. 79 80\code 81const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const 82\endcode 83 84\param[in] M a square matrix. 85\returns expression representing \f$ \cos(M) \f$. 86 87This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cos(). 88 89\sa \ref matrixbase_sin "sin()" for an example. 90 91 92 93\subsection matrixbase_cosh MatrixBase::cosh() 94 95Compute the matrix hyberbolic cosine. 96 97\code 98const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const 99\endcode 100 101\param[in] M a square matrix. 102\returns expression representing \f$ \cosh(M) \f$ 103 104This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cosh(). 105 106\sa \ref matrixbase_sinh "sinh()" for an example. 107 108 109 110\subsection matrixbase_exp MatrixBase::exp() 111 112Compute the matrix exponential. 113 114\code 115const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const 116\endcode 117 118\param[in] M matrix whose exponential is to be computed. 119\returns expression representing the matrix exponential of \p M. 120 121The matrix exponential of \f$ M \f$ is defined by 122\f[ \exp(M) = \sum_{k=0}^\infty \frac{M^k}{k!}. \f] 123The matrix exponential can be used to solve linear ordinary 124differential equations: the solution of \f$ y' = My \f$ with the 125initial condition \f$ y(0) = y_0 \f$ is given by 126\f$ y(t) = \exp(M) y_0 \f$. 127 128The cost of the computation is approximately \f$ 20 n^3 \f$ for 129matrices of size \f$ n \f$. The number 20 depends weakly on the 130norm of the matrix. 131 132The matrix exponential is computed using the scaling-and-squaring 133method combined with Padé approximation. The matrix is first 134rescaled, then the exponential of the reduced matrix is computed 135approximant, and then the rescaling is undone by repeated 136squaring. The degree of the Padé approximant is chosen such 137that the approximation error is less than the round-off 138error. However, errors may accumulate during the squaring phase. 139 140Details of the algorithm can be found in: Nicholas J. Higham, "The 141scaling and squaring method for the matrix exponential revisited," 142<em>SIAM J. %Matrix Anal. Applic.</em>, <b>26</b>:1179–1193, 1432005. 144 145Example: The following program checks that 146\f[ \exp \left[ \begin{array}{ccc} 147 0 & \frac14\pi & 0 \\ 148 -\frac14\pi & 0 & 0 \\ 149 0 & 0 & 0 150 \end{array} \right] = \left[ \begin{array}{ccc} 151 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ 152 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ 153 0 & 0 & 1 154 \end{array} \right]. \f] 155This corresponds to a rotation of \f$ \frac14\pi \f$ radians around 156the z-axis. 157 158\include MatrixExponential.cpp 159Output: \verbinclude MatrixExponential.out 160 161\note \p M has to be a matrix of \c float, \c double, \c long double 162\c complex<float>, \c complex<double>, or \c complex<long double> . 163 164 165\subsection matrixbase_log MatrixBase::log() 166 167Compute the matrix logarithm. 168 169\code 170const MatrixLogarithmReturnValue<Derived> MatrixBase<Derived>::log() const 171\endcode 172 173\param[in] M invertible matrix whose logarithm is to be computed. 174\returns expression representing the matrix logarithm root of \p M. 175 176The matrix logarithm of \f$ M \f$ is a matrix \f$ X \f$ such that 177\f$ \exp(X) = M \f$ where exp denotes the matrix exponential. As for 178the scalar logarithm, the equation \f$ \exp(X) = M \f$ may have 179multiple solutions; this function returns a matrix whose eigenvalues 180have imaginary part in the interval \f$ (-\pi,\pi] \f$. 181 182In the real case, the matrix \f$ M \f$ should be invertible and 183it should have no eigenvalues which are real and negative (pairs of 184complex conjugate eigenvalues are allowed). In the complex case, it 185only needs to be invertible. 186 187This function computes the matrix logarithm using the Schur-Parlett 188algorithm as implemented by MatrixBase::matrixFunction(). The 189logarithm of an atomic block is computed by MatrixLogarithmAtomic, 190which uses direct computation for 1-by-1 and 2-by-2 blocks and an 191inverse scaling-and-squaring algorithm for bigger blocks, with the 192square roots computed by MatrixBase::sqrt(). 193 194Details of the algorithm can be found in Section 11.6.2 of: 195Nicholas J. Higham, 196<em>Functions of Matrices: Theory and Computation</em>, 197SIAM 2008. ISBN 978-0-898716-46-7. 198 199Example: The following program checks that 200\f[ \log \left[ \begin{array}{ccc} 201 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ 202 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ 203 0 & 0 & 1 204 \end{array} \right] = \left[ \begin{array}{ccc} 205 0 & \frac14\pi & 0 \\ 206 -\frac14\pi & 0 & 0 \\ 207 0 & 0 & 0 208 \end{array} \right]. \f] 209This corresponds to a rotation of \f$ \frac14\pi \f$ radians around 210the z-axis. This is the inverse of the example used in the 211documentation of \ref matrixbase_exp "exp()". 212 213\include MatrixLogarithm.cpp 214Output: \verbinclude MatrixLogarithm.out 215 216\note \p M has to be a matrix of \c float, \c double, <tt>long 217double</tt>, \c complex<float>, \c complex<double>, or \c complex<long 218double> . 219 220\sa MatrixBase::exp(), MatrixBase::matrixFunction(), 221 class MatrixLogarithmAtomic, MatrixBase::sqrt(). 222 223 224\subsection matrixbase_pow MatrixBase::pow() 225 226Compute the matrix raised to arbitrary real power. 227 228\code 229const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(RealScalar p) const 230\endcode 231 232\param[in] M base of the matrix power, should be a square matrix. 233\param[in] p exponent of the matrix power, should be real. 234 235The matrix power \f$ M^p \f$ is defined as \f$ \exp(p \log(M)) \f$, 236where exp denotes the matrix exponential, and log denotes the matrix 237logarithm. 238 239The matrix \f$ M \f$ should meet the conditions to be an argument of 240matrix logarithm. If \p p is not of the real scalar type of \p M, it 241is casted into the real scalar type of \p M. 242 243This function computes the matrix power using the Schur-Padé 244algorithm as implemented by class MatrixPower. The exponent is split 245into integral part and fractional part, where the fractional part is 246in the interval \f$ (-1, 1) \f$. The main diagonal and the first 247super-diagonal is directly computed. 248 249Details of the algorithm can be found in: Nicholas J. Higham and 250Lijing Lin, "A Schur-Padé algorithm for fractional powers of a 251matrix," <em>SIAM J. %Matrix Anal. Applic.</em>, 252<b>32(3)</b>:1056–1078, 2011. 253 254Example: The following program checks that 255\f[ \left[ \begin{array}{ccc} 256 \cos1 & -\sin1 & 0 \\ 257 \sin1 & \cos1 & 0 \\ 258 0 & 0 & 1 259 \end{array} \right]^{\frac14\pi} = \left[ \begin{array}{ccc} 260 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ 261 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ 262 0 & 0 & 1 263 \end{array} \right]. \f] 264This corresponds to \f$ \frac14\pi \f$ rotations of 1 radian around 265the z-axis. 266 267\include MatrixPower.cpp 268Output: \verbinclude MatrixPower.out 269 270MatrixBase::pow() is user-friendly. However, there are some 271circumstances under which you should use class MatrixPower directly. 272MatrixPower can save the result of Schur decomposition, so it's 273better for computing various powers for the same matrix. 274 275Example: 276\include MatrixPower_optimal.cpp 277Output: \verbinclude MatrixPower_optimal.out 278 279\note \p M has to be a matrix of \c float, \c double, <tt>long 280double</tt>, \c complex<float>, \c complex<double>, or \c complex<long 281double> . 282 283\sa MatrixBase::exp(), MatrixBase::log(), class MatrixPower. 284 285 286\subsection matrixbase_matrixfunction MatrixBase::matrixFunction() 287 288Compute a matrix function. 289 290\code 291const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const 292\endcode 293 294\param[in] M argument of matrix function, should be a square matrix. 295\param[in] f an entire function; \c f(x,n) should compute the n-th 296derivative of f at x. 297\returns expression representing \p f applied to \p M. 298 299Suppose that \p M is a matrix whose entries have type \c Scalar. 300Then, the second argument, \p f, should be a function with prototype 301\code 302ComplexScalar f(ComplexScalar, int) 303\endcode 304where \c ComplexScalar = \c std::complex<Scalar> if \c Scalar is 305real (e.g., \c float or \c double) and \c ComplexScalar = 306\c Scalar if \c Scalar is complex. The return value of \c f(x,n) 307should be \f$ f^{(n)}(x) \f$, the n-th derivative of f at x. 308 309This routine uses the algorithm described in: 310Philip Davies and Nicholas J. Higham, 311"A Schur-Parlett algorithm for computing matrix functions", 312<em>SIAM J. %Matrix Anal. Applic.</em>, <b>25</b>:464–485, 2003. 313 314The actual work is done by the MatrixFunction class. 315 316Example: The following program checks that 317\f[ \exp \left[ \begin{array}{ccc} 318 0 & \frac14\pi & 0 \\ 319 -\frac14\pi & 0 & 0 \\ 320 0 & 0 & 0 321 \end{array} \right] = \left[ \begin{array}{ccc} 322 \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ 323 \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ 324 0 & 0 & 1 325 \end{array} \right]. \f] 326This corresponds to a rotation of \f$ \frac14\pi \f$ radians around 327the z-axis. This is the same example as used in the documentation 328of \ref matrixbase_exp "exp()". 329 330\include MatrixFunction.cpp 331Output: \verbinclude MatrixFunction.out 332 333Note that the function \c expfn is defined for complex numbers 334\c x, even though the matrix \c A is over the reals. Instead of 335\c expfn, we could also have used StdStemFunctions::exp: 336\code 337A.matrixFunction(StdStemFunctions<std::complex<double> >::exp, &B); 338\endcode 339 340 341 342\subsection matrixbase_sin MatrixBase::sin() 343 344Compute the matrix sine. 345 346\code 347const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const 348\endcode 349 350\param[in] M a square matrix. 351\returns expression representing \f$ \sin(M) \f$. 352 353This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sin(). 354 355Example: \include MatrixSine.cpp 356Output: \verbinclude MatrixSine.out 357 358 359 360\subsection matrixbase_sinh MatrixBase::sinh() 361 362Compute the matrix hyperbolic sine. 363 364\code 365MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const 366\endcode 367 368\param[in] M a square matrix. 369\returns expression representing \f$ \sinh(M) \f$ 370 371This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sinh(). 372 373Example: \include MatrixSinh.cpp 374Output: \verbinclude MatrixSinh.out 375 376 377\subsection matrixbase_sqrt MatrixBase::sqrt() 378 379Compute the matrix square root. 380 381\code 382const MatrixSquareRootReturnValue<Derived> MatrixBase<Derived>::sqrt() const 383\endcode 384 385\param[in] M invertible matrix whose square root is to be computed. 386\returns expression representing the matrix square root of \p M. 387 388The matrix square root of \f$ M \f$ is the matrix \f$ M^{1/2} \f$ 389whose square is the original matrix; so if \f$ S = M^{1/2} \f$ then 390\f$ S^2 = M \f$. 391 392In the <b>real case</b>, the matrix \f$ M \f$ should be invertible and 393it should have no eigenvalues which are real and negative (pairs of 394complex conjugate eigenvalues are allowed). In that case, the matrix 395has a square root which is also real, and this is the square root 396computed by this function. 397 398The matrix square root is computed by first reducing the matrix to 399quasi-triangular form with the real Schur decomposition. The square 400root of the quasi-triangular matrix can then be computed directly. The 401cost is approximately \f$ 25 n^3 \f$ real flops for the real Schur 402decomposition and \f$ 3\frac13 n^3 \f$ real flops for the remainder 403(though the computation time in practice is likely more than this 404indicates). 405 406Details of the algorithm can be found in: Nicholas J. Highan, 407"Computing real square roots of a real matrix", <em>Linear Algebra 408Appl.</em>, 88/89:405–430, 1987. 409 410If the matrix is <b>positive-definite symmetric</b>, then the square 411root is also positive-definite symmetric. In this case, it is best to 412use SelfAdjointEigenSolver::operatorSqrt() to compute it. 413 414In the <b>complex case</b>, the matrix \f$ M \f$ should be invertible; 415this is a restriction of the algorithm. The square root computed by 416this algorithm is the one whose eigenvalues have an argument in the 417interval \f$ (-\frac12\pi, \frac12\pi] \f$. This is the usual branch 418cut. 419 420The computation is the same as in the real case, except that the 421complex Schur decomposition is used to reduce the matrix to a 422triangular matrix. The theoretical cost is the same. Details are in: 423Åke Björck and Sven Hammarling, "A Schur method for the 424square root of a matrix", <em>Linear Algebra Appl.</em>, 42552/53:127–140, 1983. 426 427Example: The following program checks that the square root of 428\f[ \left[ \begin{array}{cc} 429 \cos(\frac13\pi) & -\sin(\frac13\pi) \\ 430 \sin(\frac13\pi) & \cos(\frac13\pi) 431 \end{array} \right], \f] 432corresponding to a rotation over 60 degrees, is a rotation over 30 degrees: 433\f[ \left[ \begin{array}{cc} 434 \cos(\frac16\pi) & -\sin(\frac16\pi) \\ 435 \sin(\frac16\pi) & \cos(\frac16\pi) 436 \end{array} \right]. \f] 437 438\include MatrixSquareRoot.cpp 439Output: \verbinclude MatrixSquareRoot.out 440 441\sa class RealSchur, class ComplexSchur, class MatrixSquareRoot, 442 SelfAdjointEigenSolver::operatorSqrt(). 443 444*/ 445 446#endif // EIGEN_MATRIX_FUNCTIONS 447 448