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