Eigenunsupported
3.3.4

This module aims to provide various methods for the computation of matrix functions.
To use this module, add
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:
Classes  
class  Eigen::MatrixComplexPowerReturnValue< Derived > 
Proxy for the matrix power of some matrix (expression). More...  
struct  Eigen::MatrixExponentialReturnValue< Derived > 
Proxy for the matrix exponential of some matrix (expression). More...  
class  Eigen::MatrixFunctionReturnValue< Derived > 
Proxy for the matrix function of some matrix (expression). More...  
class  Eigen::MatrixLogarithmReturnValue< Derived > 
Proxy for the matrix logarithm of some matrix (expression). More...  
class  Eigen::MatrixPower< MatrixType > 
Class for computing matrix powers. More...  
class  Eigen::MatrixPowerAtomic< MatrixType > 
Class for computing matrix powers. More...  
class  Eigen::MatrixPowerParenthesesReturnValue< MatrixType > 
Proxy for the matrix power of some matrix. More...  
class  Eigen::MatrixPowerReturnValue< Derived > 
Proxy for the matrix power of some matrix (expression). More...  
class  Eigen::MatrixSquareRootReturnValue< Derived > 
Proxy for the matrix square root of some matrix (expression). More...  
Functions  
template<typename MatrixType , typename ResultType >  
void  Eigen::matrix_sqrt_quasi_triangular (const MatrixType &arg, ResultType &result) 
Compute matrix square root of quasitriangular matrix. More...  
template<typename MatrixType , typename ResultType >  
void  Eigen::matrix_sqrt_triangular (const MatrixType &arg, ResultType &result) 
Compute matrix square root of triangular matrix. More...  
The remainder of the page documents the following MatrixBase methods which are defined in the MatrixFunctions module.
Compute the matrix cosine.
[in]  M  a square matrix. 
This function computes the matrix cosine. Use ArrayBase::cos() for computing the entrywise cosine.
The implementation calls matrixFunction() with StdStemFunctions::cos().
Compute the matrix hyberbolic cosine.
[in]  M  a square matrix. 
This function calls matrixFunction() with StdStemFunctions::cosh().
Compute the matrix exponential.
[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 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 for matrices of size . The number 20 depends weakly on the norm of the matrix.
The matrix exponential is computed using the scalingandsquaring 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 roundoff 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 zaxis.
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
, long
double complex<float>
, complex<double>
, or complex<long double>
.Compute the matrix logarithm.
[in]  M  invertible matrix whose logarithm is to be computed. 
M
.The matrix logarithm of is a matrix such that where exp denotes the matrix exponential. As for the scalar logarithm, the equation may have multiple solutions; this function returns a matrix whose eigenvalues have imaginary part in the interval .
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 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 SchurParlett algorithm as implemented by MatrixBase::matrixFunction(). The logarithm of an atomic block is computed by MatrixLogarithmAtomic, which uses direct computation for 1by1 and 2by2 blocks and an inverse scalingandsquaring 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, Functions of Matrices: Theory and Computation, SIAM 2008. ISBN 9780898716467.
Example: The following program checks that
This corresponds to a rotation of radians around the zaxis. This is the inverse of the example used in the documentation of exp().
Output:
The matrix A is: 0.707107 0.707107 0 0.707107 0.707107 0 0 0 1 The matrix logarithm of A is: 8.86512e17 0.785398 0 0.785398 8.86512e17 0 0 0 0
M
has to be a matrix of float
, double
, long double
, complex<float>
, complex<double>
, or complex
<long double>=""> .Compute the matrix raised to arbitrary real power.
[in]  M  base of the matrix power, should be a square matrix. 
[in]  p  exponent of the matrix power. 
The matrix power is defined as , 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 pth power. Use ArrayBase::pow() if you want to do the latter.
If p
is complex, the scalar type of M
should be the type of p
. simply evaluates into . Therefore, the matrix should meet the conditions to be an argument of matrix logarithm.
If p
is real, it is casted into the real scalar type of M
. Then this function computes the matrix power using the SchurPadé algorithm as implemented by class MatrixPower. The exponent is split into integral part and fractional part, where the fractional part is in the interval . The main diagonal and the first superdiagonal is directly computed.
If M
is singular with a semisimple zero eigenvalue and p
is positive, the Schur factor is reordered with Givens rotations, i.e.
where is invertible. Then is given by
Details of the algorithm can be found in: Nicholas J. Higham and Lijing Lin, "A SchurPadé algorithm for fractional powers of a matrix," SIAM J. Matrix Anal. Applic., 32(3):1056–1078, 2011.
Example: The following program checks that
This corresponds to rotations of 1 radian around the zaxis.
Output:
The matrix A is: 0.540302 0.841471 0 0.841471 0.540302 0 0 0 1 The matrix power A^(pi/4) is: 0.707107 0.707107 0 0.707107 0.707107 0 0 0 1
MatrixBase::pow() is userfriendly. 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:
Output:
The matrix A is: (0.211234,0.680375) (0.10794,0.444451) (0.434594,0.271423) (0.198111,0.686642) (0.59688,0.566198) (0.257742,0.0452059) (0.213938,0.716795) (0.782382,0.740419) (0.604897,0.823295) (0.0268018,0.270431) (0.514226,0.967399) (0.563486,0.997849) (0.536459,0.329554) (0.83239,0.904459) (0.608354,0.725537) (0.678224,0.0258648) A^3.1 is: (2.80575,0.607662) (1.16847,0.00660555) (0.760385,1.01461) (0.38073,0.106512) (1.4041,3.61891) (1.00481,0.186263) (0.163888,0.449419) (0.388981,1.22629) (2.07957,1.58136) (0.825866,2.25962) (5.09383,0.155736) (0.394308,1.63034) (0.818997,0.671026) (2.11069,0.00768024) (1.37876,0.140165) (2.50512,0.854429) A^3.3 is: (2.83571,0.238717) (1.48174,0.0615217) (0.0544396,1.68092) (0.292699,0.621726) (2.0521,3.58316) (0.87894,0.400548) (0.738072,0.121242) (1.07957,1.63492) (3.00106,1.10558) (1.52205,1.92407) (5.29759,1.83562) (0.532038,1.50253) (0.491353,0.4145) (2.5761,0.481286) (1.21994,0.0367069) (2.67112,1.06331) A^3.7 is: (1.42126,0.33362) (1.39486,0.560486) (1.44968,2.47066) (0.324079,1.75879) (2.65301,1.82427) (0.357333,0.192429) (2.01017,1.4791) (2.71518,2.35892) (3.98544,0.964861) (2.26033,0.554254) (3.18211,5.94352) (2.22888,0.128951) (0.944969,2.14683) (3.31345,1.66075) (0.0623743,0.848324) (2.3897,1.863) A^3.9 is: (0.0720766,0.378685) (0.931961,0.978624) (1.9855,2.34105) (0.530547,2.17664) (2.40934,0.265286) (0.0299975,1.08827) (1.98974,2.05886) (3.45767,2.50235) (3.71666,2.3874) (2.054,0.303) (0.844348,7.29588) (2.59136,1.57689) (1.87645,2.38798) (3.52111,2.10508) (0.799055,1.6122) (1.93452,2.44408)
M
has to be a matrix of float
, double
, long double
, complex<float>
, complex<double>
, or complex
<long double>=""> .Compute a matrix function.
[in]  M  argument of matrix function, should be a square matrix. 
[in]  f  an entire function; f(x,n) should compute the nth 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
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 nth derivative of f at x.
This routine uses the algorithm described in: Philip Davies and Nicholas J. Higham, "A SchurParlett 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 zaxis. This is the same example as used in the documentation of exp().
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:
Compute the matrix sine.
[in]  M  a square matrix. 
This function computes the matrix sine. Use ArrayBase::sin() for computing the entrywise sine.
The implementation calls matrixFunction() with StdStemFunctions::sin().
Example:
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 7.77156e16 4.71845e16 5.55112e17 1 2.77556e16 1.66533e16 2.08167e16 1
Compute the matrix hyperbolic sine.
[in]  M  a square matrix. 
This function calls matrixFunction() with StdStemFunctions::sinh().
Example:
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.739988 0.256871 0.194928 0.826512 0.537546 0.562584 0.53163 0.425199 cosh(A) = 1.07817 0.567068 0.132125 0.004186 1.11649 0.135361 0.128891 0.065999 0.851201 cosh^2(A)  sinh^2(A) = 1 1.19209e07 0 2.10479e07 1 1.78814e07 1.19209e07 2.38419e07 1
Compute the matrix square root.
[in]  M  invertible matrix whose square root is to be computed. 
M
.The matrix square root of is the matrix whose square is the original matrix; so if then . 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 real case, the matrix 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 quasitriangular form with the real Schur decomposition. The square root of the quasitriangular matrix can then be computed directly. The cost is approximately real flops for the real Schur decomposition and 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", Linear Algebra Appl., 88/89:405–430, 1987.
If the matrix is positivedefinite symmetric, then the square root is also positivedefinite symmetric. In this case, it is best to use SelfAdjointEigenSolver::operatorSqrt() to compute it.
In the complex case, the matrix 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 . 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", Linear Algebra Appl., 52/53:127–140, 1983.
Example: The following program checks that the square root of
corresponding to a rotation over 60 degrees, is a rotation over 30 degrees:
Output:
The matrix A is: 0.5 0.866025 0.866025 0.5 The matrix square root of A is: 0.866025 0.5 0.5 0.866025 The square of the last matrix is: 0.5 0.866025 0.866025 0.5
void Eigen::matrix_sqrt_quasi_triangular  (  const MatrixType &  arg, 
ResultType &  result  
) 
Compute matrix square root of quasitriangular matrix.
MatrixType  type of arg , the argument of matrix square root, expected to be an instantiation of the Matrix class template. 
ResultType  type of result , where result is to be stored. 
[in]  arg  argument of matrix square root. 
[out]  result  matrix square root of upper Hessenberg part of arg . 
This function computes the square root of the upper quasitriangular matrix stored in the upper Hessenberg part of arg
. Only the upper Hessenberg part of result
is updated, the rest is not touched. See MatrixBase::sqrt() for details on how this computation is implemented.
void Eigen::matrix_sqrt_triangular  (  const MatrixType &  arg, 
ResultType &  result  
) 
Compute matrix square root of triangular matrix.
MatrixType  type of arg , the argument of matrix square root, expected to be an instantiation of the Matrix class template. 
ResultType  type of result , where result is to be stored. 
[in]  arg  argument of matrix square root. 
[out]  result  matrix square root of upper triangular part of arg . 
Only the upper triangular part (including the diagonal) of result
is updated, the rest is not touched. See MatrixBase::sqrt() for details on how this computation is implemented.