501 lines
17 KiB
Plaintext
501 lines
17 KiB
Plaintext
// This file is part of Eigen, a lightweight C++ template library
|
|
// for linear algebra.
|
|
//
|
|
// Copyright (C) 2009 Jitse Niesen <jitse@maths.leeds.ac.uk>
|
|
// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
|
|
//
|
|
// This Source Code Form is subject to the terms of the Mozilla
|
|
// Public License v. 2.0. If a copy of the MPL was not distributed
|
|
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
|
|
|
#ifndef EIGEN_MATRIX_FUNCTIONS
|
|
#define EIGEN_MATRIX_FUNCTIONS
|
|
|
|
#include <cfloat>
|
|
#include <list>
|
|
|
|
#include <Eigen/Core>
|
|
#include <Eigen/LU>
|
|
#include <Eigen/Eigenvalues>
|
|
|
|
/**
|
|
* \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_log "MatrixBase::log()", for computing the matrix logarithm
|
|
* - \ref matrixbase_pow "MatrixBase::pow()", for computing the matrix power
|
|
* - \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"
|
|
#include "src/MatrixFunctions/MatrixLogarithm.h"
|
|
#include "src/MatrixFunctions/MatrixPower.h"
|
|
|
|
|
|
/**
|
|
\page matrixbaseextra_page
|
|
\ingroup MatrixFunctions_Module
|
|
|
|
\section matrixbaseextra MatrixBase methods defined in the MatrixFunctions module
|
|
|
|
The remainder of the page documents the following MatrixBase methods
|
|
which are defined in the MatrixFunctions module.
|
|
|
|
|
|
|
|
\subsection 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 computes the matrix cosine. Use ArrayBase::cos() for computing the entry-wise cosine.
|
|
|
|
The implementation calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cos().
|
|
|
|
\sa \ref matrixbase_sin "sin()" for an example.
|
|
|
|
|
|
|
|
\subsection 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.
|
|
|
|
|
|
|
|
\subsection 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 matrix exponential is different from applying the exp function to all the entries in the matrix.
|
|
Use ArrayBase::exp() if you want to do the latter.
|
|
|
|
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é 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é 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–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, `long double`
|
|
\c complex<float>, \c complex<double>, or `complex<long double>` .
|
|
|
|
|
|
\subsection matrixbase_log MatrixBase::log()
|
|
|
|
Compute the matrix logarithm.
|
|
|
|
\code
|
|
const MatrixLogarithmReturnValue<Derived> MatrixBase<Derived>::log() const
|
|
\endcode
|
|
|
|
\param[in] M invertible matrix whose logarithm is to be computed.
|
|
\returns expression representing the matrix logarithm root of \p M.
|
|
|
|
The matrix logarithm of \f$ M \f$ is a matrix \f$ X \f$ such that
|
|
\f$ \exp(X) = M \f$ where exp denotes the matrix exponential. As for
|
|
the scalar logarithm, the equation \f$ \exp(X) = M \f$ may have
|
|
multiple solutions; this function returns a matrix whose eigenvalues
|
|
have imaginary part in the interval \f$ (-\pi,\pi] \f$.
|
|
|
|
The matrix logarithm is different from applying the log function to all the entries in the matrix.
|
|
Use ArrayBase::log() if you want to do the latter.
|
|
|
|
In the real case, 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 the complex case, it
|
|
only needs to be invertible.
|
|
|
|
This function computes the matrix logarithm using the Schur-Parlett
|
|
algorithm as implemented by MatrixBase::matrixFunction(). The
|
|
logarithm of an atomic block is computed by MatrixLogarithmAtomic,
|
|
which uses direct computation for 1-by-1 and 2-by-2 blocks and an
|
|
inverse scaling-and-squaring algorithm for bigger blocks, with the
|
|
square roots computed by MatrixBase::sqrt().
|
|
|
|
Details of the algorithm can be found in Section 11.6.2 of:
|
|
Nicholas J. Higham,
|
|
<em>Functions of Matrices: Theory and Computation</em>,
|
|
SIAM 2008. ISBN 978-0-898716-46-7.
|
|
|
|
Example: The following program checks that
|
|
\f[ \log \left[ \begin{array}{ccc}
|
|
\frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
|
|
\frac12\sqrt2 & \frac12\sqrt2 & 0 \\
|
|
0 & 0 & 1
|
|
\end{array} \right] = \left[ \begin{array}{ccc}
|
|
0 & \frac14\pi & 0 \\
|
|
-\frac14\pi & 0 & 0 \\
|
|
0 & 0 & 0
|
|
\end{array} \right]. \f]
|
|
This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
|
|
the z-axis. This is the inverse of the example used in the
|
|
documentation of \ref matrixbase_exp "exp()".
|
|
|
|
\include MatrixLogarithm.cpp
|
|
Output: \verbinclude MatrixLogarithm.out
|
|
|
|
\note \p M has to be a matrix of \c float, \c double, `long
|
|
double`, \c complex<float>, \c complex<double>, or `complex<long double>`.
|
|
|
|
\sa MatrixBase::exp(), MatrixBase::matrixFunction(),
|
|
class MatrixLogarithmAtomic, MatrixBase::sqrt().
|
|
|
|
|
|
\subsection matrixbase_pow MatrixBase::pow()
|
|
|
|
Compute the matrix raised to arbitrary real power.
|
|
|
|
\code
|
|
const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(RealScalar p) const
|
|
\endcode
|
|
|
|
\param[in] M base of the matrix power, should be a square matrix.
|
|
\param[in] p exponent of the matrix power.
|
|
|
|
The matrix power \f$ M^p \f$ is defined as \f$ \exp(p \log(M)) \f$,
|
|
where exp denotes the matrix exponential, and log denotes the matrix
|
|
logarithm. This is different from raising all the entries in the matrix
|
|
to the p-th power. Use ArrayBase::pow() if you want to do the latter.
|
|
|
|
If \p p is complex, the scalar type of \p M should be the type of \p
|
|
p . \f$ M^p \f$ simply evaluates into \f$ \exp(p \log(M)) \f$.
|
|
Therefore, the matrix \f$ M \f$ should meet the conditions to be an
|
|
argument of matrix logarithm.
|
|
|
|
If \p p is real, it is casted into the real scalar type of \p M. Then
|
|
this function computes the matrix power using the Schur-Padé
|
|
algorithm as implemented by class MatrixPower. The exponent is split
|
|
into integral part and fractional part, where the fractional part is
|
|
in the interval \f$ (-1, 1) \f$. The main diagonal and the first
|
|
super-diagonal is directly computed.
|
|
|
|
If \p M is singular with a semisimple zero eigenvalue and \p p is
|
|
positive, the Schur factor \f$ T \f$ is reordered with Givens
|
|
rotations, i.e.
|
|
|
|
\f[ T = \left[ \begin{array}{cc}
|
|
T_1 & T_2 \\
|
|
0 & 0
|
|
\end{array} \right] \f]
|
|
|
|
where \f$ T_1 \f$ is invertible. Then \f$ T^p \f$ is given by
|
|
|
|
\f[ T^p = \left[ \begin{array}{cc}
|
|
T_1^p & T_1^{-1} T_1^p T_2 \\
|
|
0 & 0
|
|
\end{array}. \right] \f]
|
|
|
|
\warning Fractional power of a matrix with a non-semisimple zero
|
|
eigenvalue is not well-defined. We introduce an assertion failure
|
|
against inaccurate result, e.g. \code
|
|
#include <unsupported/Eigen/MatrixFunctions>
|
|
#include <iostream>
|
|
|
|
int main()
|
|
{
|
|
Eigen::Matrix4d A;
|
|
A << 0, 0, 2, 3,
|
|
0, 0, 4, 5,
|
|
0, 0, 6, 7,
|
|
0, 0, 8, 9;
|
|
std::cout << A.pow(0.37) << std::endl;
|
|
|
|
// The 1 makes eigenvalue 0 non-semisimple.
|
|
A.coeffRef(0, 1) = 1;
|
|
|
|
// This fails if EIGEN_NO_DEBUG is undefined.
|
|
std::cout << A.pow(0.37) << std::endl;
|
|
|
|
return 0;
|
|
}
|
|
\endcode
|
|
|
|
Details of the algorithm can be found in: Nicholas J. Higham and
|
|
Lijing Lin, "A Schur-Padé algorithm for fractional powers of a
|
|
matrix," <em>SIAM J. %Matrix Anal. Applic.</em>,
|
|
<b>32(3)</b>:1056–1078, 2011.
|
|
|
|
Example: The following program checks that
|
|
\f[ \left[ \begin{array}{ccc}
|
|
\cos1 & -\sin1 & 0 \\
|
|
\sin1 & \cos1 & 0 \\
|
|
0 & 0 & 1
|
|
\end{array} \right]^{\frac14\pi} = \left[ \begin{array}{ccc}
|
|
\frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
|
|
\frac12\sqrt2 & \frac12\sqrt2 & 0 \\
|
|
0 & 0 & 1
|
|
\end{array} \right]. \f]
|
|
This corresponds to \f$ \frac14\pi \f$ rotations of 1 radian around
|
|
the z-axis.
|
|
|
|
\include MatrixPower.cpp
|
|
Output: \verbinclude MatrixPower.out
|
|
|
|
MatrixBase::pow() is user-friendly. However, there are some
|
|
circumstances under which you should use class MatrixPower directly.
|
|
MatrixPower can save the result of Schur decomposition, so it's
|
|
better for computing various powers for the same matrix.
|
|
|
|
Example:
|
|
\include MatrixPower_optimal.cpp
|
|
Output: \verbinclude MatrixPower_optimal.out
|
|
|
|
\note \p M has to be a matrix of \c float, \c double, `long
|
|
double`, \c complex<float>, \c complex<double>, or
|
|
\c complex<long double> .
|
|
|
|
\sa MatrixBase::exp(), MatrixBase::log(), class MatrixPower.
|
|
|
|
|
|
\subsection 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–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
|
|
|
|
|
|
|
|
\subsection 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 computes the matrix sine. Use ArrayBase::sin() for computing the entry-wise sine.
|
|
|
|
The implementation calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sin().
|
|
|
|
Example: \include MatrixSine.cpp
|
|
Output: \verbinclude MatrixSine.out
|
|
|
|
|
|
|
|
\subsection 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
|
|
|
|
|
|
\subsection 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$. This is different from taking the square root of all
|
|
the entries in the matrix; use ArrayBase::sqrt() if you want to do the
|
|
latter.
|
|
|
|
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–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:
|
|
Åke Björck and Sven Hammarling, "A Schur method for the
|
|
square root of a matrix", <em>Linear Algebra Appl.</em>,
|
|
52/53:127–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
|
|
|