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