eigen/unsupported/Eigen/MatrixFunctions

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.

#ifndef EIGEN_MATRIX_FUNCTIONS
#define EIGEN_MATRIX_FUNCTIONS

#include <list>
#include <functional>
#include <iterator>

#include <Eigen/Core>
#include <Eigen/LU>
#include <Eigen/Eigenvalues>

namespace Eigen {

/** \ingroup Unsupported_modules
  * \defgroup MatrixFunctions_Module Matrix functions module
  * \brief This module aims to provide various methods for the computation of
  * matrix functions. 
  *
  * To use this module, add 
  * \code
  * #include <unsupported/Eigen/MatrixFunctions>
  * \endcode
  * at the start of your source file.
  *
  * This module defines the following MatrixBase methods.
  *  - \ref matrixbase_cos "MatrixBase::cos()", for computing the matrix cosine
  *  - \ref matrixbase_cosh "MatrixBase::cosh()", for computing the matrix hyperbolic cosine
  *  - \ref matrixbase_exp "MatrixBase::exp()", for computing the matrix exponential
  *  - \ref matrixbase_matrixfunction "MatrixBase::matrixFunction()", for computing general matrix functions
  *  - \ref matrixbase_sin "MatrixBase::sin()", for computing the matrix sine
  *  - \ref matrixbase_sinh "MatrixBase::sinh()", for computing the matrix hyperbolic sine
  *  - \ref matrixbase_sqrt "MatrixBase::sqrt()", for computing the matrix square root
  *
  * These methods are the main entry points to this module. 
  *
  * %Matrix functions are defined as follows.  Suppose that \f$ f \f$
  * is an entire function (that is, a function on the complex plane
  * that is everywhere complex differentiable).  Then its Taylor
  * series
  * \f[ f(0) + f'(0) x + \frac{f''(0)}{2} x^2 + \frac{f'''(0)}{3!} x^3 + \cdots \f]
  * converges to \f$ f(x) \f$. In this case, we can define the matrix
  * function by the same series:
  * \f[ f(M) = f(0) + f'(0) M + \frac{f''(0)}{2} M^2 + \frac{f'''(0)}{3!} M^3 + \cdots \f]
  *
  */

#include "src/MatrixFunctions/MatrixExponential.h"
#include "src/MatrixFunctions/MatrixFunction.h"
#include "src/MatrixFunctions/MatrixSquareRoot.h"



/** 
\page matrixbaseextra MatrixBase methods defined in the MatrixFunctions module
\ingroup MatrixFunctions_Module

The remainder of the page documents the following MatrixBase methods
which are defined in the MatrixFunctions module.



\section matrixbase_cos MatrixBase::cos()

Compute the matrix cosine.

\code
const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
\endcode

\param[in]  M  a square matrix.
\returns  expression representing \f$ \cos(M) \f$.

This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cos().

\sa \ref matrixbase_sin "sin()" for an example.



\section matrixbase_cosh MatrixBase::cosh()

Compute the matrix hyberbolic cosine.

\code
const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const
\endcode

\param[in]  M  a square matrix.
\returns  expression representing \f$ \cosh(M) \f$

This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cosh().

\sa \ref matrixbase_sinh "sinh()" for an example.



\section matrixbase_exp MatrixBase::exp()

Compute the matrix exponential.

\code
const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
\endcode

\param[in]  M  matrix whose exponential is to be computed.
\returns    expression representing the matrix exponential of \p M.

The matrix exponential of \f$ M \f$ is defined by
\f[ \exp(M) = \sum_{k=0}^\infty \frac{M^k}{k!}. \f]
The matrix exponential can be used to solve linear ordinary
differential equations: the solution of \f$ y' = My \f$ with the
initial condition \f$ y(0) = y_0 \f$ is given by
\f$ y(t) = \exp(M) y_0 \f$.

The cost of the computation is approximately \f$ 20 n^3 \f$ for
matrices of size \f$ n \f$. The number 20 depends weakly on the
norm of the matrix.

The matrix exponential is computed using the scaling-and-squaring
method combined with Pad&eacute; approximation. The matrix is first
rescaled, then the exponential of the reduced matrix is computed
approximant, and then the rescaling is undone by repeated
squaring. The degree of the Pad&eacute; approximant is chosen such
that the approximation error is less than the round-off
error. However, errors may accumulate during the squaring phase.

Details of the algorithm can be found in: Nicholas J. Higham, "The
scaling and squaring method for the matrix exponential revisited,"
<em>SIAM J. %Matrix Anal. Applic.</em>, <b>26</b>:1179&ndash;1193,
2005.

Example: The following program checks that
\f[ \exp \left[ \begin{array}{ccc}
      0 & \frac14\pi & 0 \\
      -\frac14\pi & 0 & 0 \\
      0 & 0 & 0
    \end{array} \right] = \left[ \begin{array}{ccc}
      \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
      \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
      0 & 0 & 1
    \end{array} \right]. \f]
This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
the z-axis.

\include MatrixExponential.cpp
Output: \verbinclude MatrixExponential.out

\note \p M has to be a matrix of \c float, \c double,
\c complex<float> or \c complex<double> .



\section matrixbase_matrixfunction MatrixBase::matrixFunction()

Compute a matrix function.

\code
const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const
\endcode

\param[in]  M  argument of matrix function, should be a square matrix.
\param[in]  f  an entire function; \c f(x,n) should compute the n-th
derivative of f at x.
\returns  expression representing \p f applied to \p M.

Suppose that \p M is a matrix whose entries have type \c Scalar. 
Then, the second argument, \p f, should be a function with prototype
\code 
ComplexScalar f(ComplexScalar, int) 
\endcode
where \c ComplexScalar = \c std::complex<Scalar> if \c Scalar is
real (e.g., \c float or \c double) and \c ComplexScalar =
\c Scalar if \c Scalar is complex. The return value of \c f(x,n)
should be \f$ f^{(n)}(x) \f$, the n-th derivative of f at x.

This routine uses the algorithm described in:
Philip Davies and Nicholas J. Higham, 
"A Schur-Parlett algorithm for computing matrix functions", 
<em>SIAM J. %Matrix Anal. Applic.</em>, <b>25</b>:464&ndash;485, 2003.

The actual work is done by the MatrixFunction class.

Example: The following program checks that
\f[ \exp \left[ \begin{array}{ccc} 
      0 & \frac14\pi & 0 \\ 
      -\frac14\pi & 0 & 0 \\
      0 & 0 & 0 
    \end{array} \right] = \left[ \begin{array}{ccc}
      \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
      \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
      0 & 0 & 1
    \end{array} \right]. \f]
This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
the z-axis. This is the same example as used in the documentation
of \ref matrixbase_exp "exp()".

\include MatrixFunction.cpp
Output: \verbinclude MatrixFunction.out

Note that the function \c expfn is defined for complex numbers 
\c x, even though the matrix \c A is over the reals. Instead of
\c expfn, we could also have used StdStemFunctions::exp:
\code
A.matrixFunction(StdStemFunctions<std::complex<double> >::exp, &B);
\endcode



\section matrixbase_sin MatrixBase::sin()

Compute the matrix sine.

\code
const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
\endcode

\param[in]  M  a square matrix.
\returns  expression representing \f$ \sin(M) \f$.

This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sin().

Example: \include MatrixSine.cpp
Output: \verbinclude MatrixSine.out



\section matrixbase_sinh MatrixBase::sinh()

Compute the matrix hyperbolic sine.

\code
MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const
\endcode

\param[in]  M  a square matrix.
\returns  expression representing \f$ \sinh(M) \f$

This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sinh().

Example: \include MatrixSinh.cpp
Output: \verbinclude MatrixSinh.out


\section matrixbase_sqrt MatrixBase::sqrt()

Compute the matrix square root.

\code
const MatrixSquareRootReturnValue<Derived> MatrixBase<Derived>::sqrt() const
\endcode

\param[in]  M  invertible matrix whose square root is to be computed.
\returns    expression representing the matrix square root of \p M.

The matrix square root of \f$ M \f$ is the matrix \f$ M^{1/2} \f$
whose square is the original matrix; so if \f$ S = M^{1/2} \f$ then
\f$ S^2 = M \f$. 

In the <b>real case</b>, the matrix \f$ M \f$ should be invertible and
it should have no eigenvalues which are real and negative (pairs of
complex conjugate eigenvalues are allowed). In that case, the matrix
has a square root which is also real, and this is the square root
computed by this function. 

The matrix square root is computed by first reducing the matrix to
quasi-triangular form with the real Schur decomposition. The square
root of the quasi-triangular matrix can then be computed directly. The
cost is approximately \f$ 25 n^3 \f$ real flops for the real Schur
decomposition and \f$ 3\frac13 n^3 \f$ real flops for the remainder
(though the computation time in practice is likely more than this
indicates).

Details of the algorithm can be found in: Nicholas J. Highan,
"Computing real square roots of a real matrix", <em>Linear Algebra
Appl.</em>, 88/89:405&ndash;430, 1987.

If the matrix is <b>positive-definite symmetric</b>, then the square
root is also positive-definite symmetric. In this case, it is best to
use SelfAdjointEigenSolver::operatorSqrt() to compute it.

In the <b>complex case</b>, the matrix \f$ M \f$ should be invertible;
this is a restriction of the algorithm. The square root computed by
this algorithm is the one whose eigenvalues have an argument in the
interval \f$ (-\frac12\pi, \frac12\pi] \f$. This is the usual branch
cut.

The computation is the same as in the real case, except that the
complex Schur decomposition is used to reduce the matrix to a
triangular matrix. The theoretical cost is the same. Details are in:
&Aring;ke Bj&ouml;rck and Sven Hammarling, "A Schur method for the
square root of a matrix", <em>Linear Algebra Appl.</em>,
52/53:127&ndash;140, 1983.

Example: The following program checks that the square root of
\f[ \left[ \begin{array}{cc} 
              \cos(\frac13\pi) & -\sin(\frac13\pi) \\
              \sin(\frac13\pi) & \cos(\frac13\pi)
    \end{array} \right], \f]
corresponding to a rotation over 60 degrees, is a rotation over 30 degrees:
\f[ \left[ \begin{array}{cc} 
              \cos(\frac16\pi) & -\sin(\frac16\pi) \\
              \sin(\frac16\pi) & \cos(\frac16\pi)
    \end{array} \right]. \f]

\include MatrixSquareRoot.cpp
Output: \verbinclude MatrixSquareRoot.out

\sa class RealSchur, class ComplexSchur, class MatrixSquareRoot,
    SelfAdjointEigenSolver::operatorSqrt().

*/

}

#endif // EIGEN_MATRIX_FUNCTIONS