This module aims to provide various methods for the computation of matrix functions. More...
Classes | |
class | MatrixExponential< MatrixType > |
Class for computing the matrix exponential. More... | |
struct | MatrixExponentialReturnValue< Derived > |
Proxy for the matrix exponential of some matrix (expression). More... | |
class | MatrixFunction< MatrixType, IsComplex > |
Class for computing matrix exponentials. More... | |
class | MatrixFunction< MatrixType, 0 > |
Partial specialization of MatrixFunction for real matrices. More... | |
class | MatrixFunction< MatrixType, 1 > |
Partial specialization of MatrixFunction for complex matrices. More... | |
class | MatrixFunctionAtomic< MatrixType > |
Helper class for computing matrix functions of atomic matrices. More... | |
class | MatrixFunctionReturnValue< Derived > |
Proxy for the matrix function of some matrix (expression). More... | |
class | StdStemFunctions< Scalar > |
Stem functions corresponding to standard mathematical functions. More... |
This module aims to provide various methods for the computation of matrix functions.
To use this module, add
#include <unsupported/Eigen/MatrixFunctions>
at the start of your source file.
This module defines the following MatrixBase methods.
These methods are the main entry points to this module.
Matrix functions are defined as follows. Suppose that is an entire function (that is, a function on the complex plane that is everywhere complex differentiable). Then its Taylor series
converges to . In this case, we can define the matrix function by the same series:
The remainder of the page documents the following MatrixBase methods which are defined in the MatrixFunctions module.
Compute the matrix cosine.
const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
[in] | M | a square matrix. |
This function calls matrixFunction() with StdStemFunctions::cos().
Compute the matrix hyberbolic cosine.
const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const
[in] | M | a square matrix. |
This function calls matrixFunction() with StdStemFunctions::cosh().
Compute the matrix exponential.
const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
[in] | M | matrix whose exponential is to be computed. |
M
.The matrix exponential of is defined by
The matrix exponential can be used to solve linear ordinary differential equations: the solution of with the initial condition
is given by
.
The cost of the computation is approximately for matrices of size
. 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," SIAM J. Matrix Anal. Applic., 26:1179–1193, 2005.
Example: The following program checks that
This corresponds to a rotation of radians around the z-axis.
#include <unsupported/Eigen/MatrixFunctions> #include <iostream> using namespace Eigen; int main() { const double pi = std::acos(-1.0); MatrixXd A(3,3); A << 0, -pi/4, 0, pi/4, 0, 0, 0, 0, 0; std::cout << "The matrix A is:\n" << A << "\n\n"; std::cout << "The matrix exponential of A is:\n" << A.exp() << "\n\n"; }
Output:
The matrix A is: 0 -0.785398 0 0.785398 0 0 0 0 0 The matrix exponential of A is: 0.707107 -0.707107 0 0.707107 0.707107 0 0 0 1
M
has to be a matrix of float
, double
, complex<float>
or complex<double>
.Compute a matrix function.
const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const
[in] | M | argument of matrix function, should be a square matrix. |
[in] | f | an entire function; f(x,n) should compute the n-th derivative of f at x. |
f
applied to M
.Suppose that M
is a matrix whose entries have type Scalar
. Then, the second argument, f
, should be a function with prototype
ComplexScalar f(ComplexScalar, int)
where ComplexScalar
= std::complex<Scalar>
if Scalar
is real (e.g., float
or double
) and ComplexScalar
= Scalar
if Scalar
is complex. The return value of f(x,n)
should be , 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", SIAM J. Matrix Anal. Applic., 25:464–485, 2003.
The actual work is done by the MatrixFunction class.
Example: The following program checks that
This corresponds to a rotation of radians around the z-axis. This is the same example as used in the documentation of exp().
#include <unsupported/Eigen/MatrixFunctions> #include <iostream> using namespace Eigen; std::complex<double> expfn(std::complex<double> x, int) { return std::exp(x); } int main() { const double pi = std::acos(-1.0); MatrixXd A(3,3); A << 0, -pi/4, 0, pi/4, 0, 0, 0, 0, 0; std::cout << "The matrix A is:\n" << A << "\n\n"; std::cout << "The matrix exponential of A is:\n" << A.matrixFunction(expfn) << "\n\n"; }
Output:
The matrix A is: 0 -0.785398 0 0.785398 0 0 0 0 0 The matrix exponential of A is: 0.707107 -0.707107 0 0.707107 0.707107 0 0 0 1
Note that the function expfn
is defined for complex numbers x
, even though the matrix A
is over the reals. Instead of expfn
, we could also have used StdStemFunctions::exp:
A.matrixFunction(StdStemFunctions<std::complex<double> >::exp, &B);
Compute the matrix sine.
const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
[in] | M | a square matrix. |
This function calls matrixFunction() with StdStemFunctions::sin().
Example:
#include <unsupported/Eigen/MatrixFunctions> #include <iostream> using namespace Eigen; int main() { MatrixXd A = MatrixXd::Random(3,3); std::cout << "A = \n" << A << "\n\n"; MatrixXd sinA = A.sin(); std::cout << "sin(A) = \n" << sinA << "\n\n"; MatrixXd cosA = A.cos(); std::cout << "cos(A) = \n" << cosA << "\n\n"; // The matrix functions satisfy sin^2(A) + cos^2(A) = I, // like the scalar functions. std::cout << "sin^2(A) + cos^2(A) = \n" << sinA*sinA + cosA*cosA << "\n\n"; }
Output:
A = 0.680375 0.59688 -0.329554 -0.211234 0.823295 0.536459 0.566198 -0.604897 -0.444451 sin(A) = 0.679919 0.4579 -0.400612 -0.227278 0.821913 0.5358 0.570141 -0.676728 -0.462398 cos(A) = 0.927728 -0.530361 -0.110482 0.00969246 0.889022 -0.137604 -0.132574 -0.04289 1.16475 sin^2(A) + cos^2(A) = 1 2.22045e-16 4.44089e-16 4.02456e-16 1 -5.55112e-16 -3.33067e-16 -3.57353e-16 1
Compute the matrix hyperbolic sine.
MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const
[in] | M | a square matrix. |
This function calls matrixFunction() with StdStemFunctions::sinh().
Example:
#include <unsupported/Eigen/MatrixFunctions> #include <iostream> using namespace Eigen; int main() { MatrixXf A = MatrixXf::Random(3,3); std::cout << "A = \n" << A << "\n\n"; MatrixXf sinhA = A.sinh(); std::cout << "sinh(A) = \n" << sinhA << "\n\n"; MatrixXf coshA = A.cosh(); std::cout << "cosh(A) = \n" << coshA << "\n\n"; // The matrix functions satisfy cosh^2(A) - sinh^2(A) = I, // like the scalar functions. std::cout << "cosh^2(A) - sinh^2(A) = \n" << coshA*coshA - sinhA*sinhA << "\n\n"; }
Output:
A = 0.680375 0.59688 -0.329554 -0.211234 0.823295 0.536459 0.566198 -0.604897 -0.444451 sinh(A) = 0.682534 0.739989 -0.256871 -0.194928 0.826512 0.537546 0.562585 -0.53163 -0.425199 cosh(A) = 1.07817 0.567068 0.132125 -0.00418616 1.11649 0.135361 0.128891 0.0659992 0.851201 cosh^2(A) - sinh^2(A) = 1 1.19209e-07 1.78814e-07 -5.58794e-08 1 0 1.19209e-07 -1.19209e-07 1