Spec-Zone.ru › Eigen3

Модуль функций матриц

Этот модуль предназначен для предоставления различных методов вычисления функций матриц.

Для использования этого модуля добавьте

#include <unsupported/Eigen/MatrixFunctions>

в начало вашего исходного файла.

Этот модуль определяет следующие методы MatrixBase.

  • MatrixBase::cos(), для вычисления косинуса матрицы
  • MatrixBase::cosh(), для вычисления гиперболического косинуса матрицы
  • MatrixBase::exp(), для вычисления экспоненты матрицы
  • MatrixBase::log(), для вычисления логарифма матрицы
  • MatrixBase::pow(), для вычисления степени матрицы
  • MatrixBase::matrixFunction(), для вычисления общих функций матриц
  • MatrixBase::sin(), для вычисления синуса матрицы
  • MatrixBase::sinh(), для вычисления гиперболического синуса матрицы
  • MatrixBase::sqrt(), для вычисления квадратного корня из матрицы

Эти методы являются основными точками входа в этот модуль.

Функции матриц определяются следующим образом. Предположим, что \( f \) — целая функция (то есть функция на комплексной плоскости, которая везде комплексно дифференцируема). Тогда её ряд Тейлора

\[ f(0) + f'(0) x + \frac{f''(0)}{2} x^2 + \frac{f'''(0)}{3!} x^3 + \cdots \]

сходится к \( f(x) \). В этом случае мы можем определить функцию матрицы тем же рядом:

\[ f(M) = f(0) + f'(0) M + \frac{f''(0)}{2} M^2 + \frac{f'''(0)}{3!} M^3 + \cdots \]

класс Eigen::MatrixComplexPowerReturnValue< Derived >
Прокси для степени матрицы некоторой матрицы (выражения). Подробнее...
структура Eigen::MatrixExponentialReturnValue< Derived >
Прокси для экспоненты матрицы некоторой матрицы (выражения). Подробнее...
класс Eigen::MatrixFunctionReturnValue< Derived >
Прокси для функции матрицы некоторой матрицы (выражения). Подробнее...
класс Eigen::MatrixLogarithmReturnValue< Derived >
Прокси для логарифма матрицы некоторой матрицы (выражения). Подробнее...
класс Eigen::MatrixPower< MatrixType >
Класс для вычисления степеней матриц. Подробнее...
класс Eigen::MatrixPowerAtomic< MatrixType >
Класс для вычисления степеней матриц. Подробнее...
класс Eigen::MatrixPowerParenthesesReturnValue< MatrixType >
Прокси для степени матрицы некоторой матрицы. Подробнее...
класс Eigen::MatrixPowerReturnValue< Derived >
Прокси для степени матрицы некоторой матрицы (выражения). Подробнее...
класс Eigen::MatrixSquareRootReturnValue< Derived >
Прокси для квадратного корня из матрицы некоторой матрицы (выражения). Подробнее...
шаблон<typename MatrixType , typename ResultType >
void Eigen::matrix_sqrt_quasi_triangular (const MatrixType &arg, ResultType &result)
Вычисление квадратного корня из квазитреугольной матрицы. Подробнее...
шаблон<typename MatrixType , typename ResultType >
void Eigen::matrix_sqrt_triangular (const MatrixType &arg, ResultType &result)
Вычисление квадратного корня из треугольной матрицы. Подробнее...

Методы MatrixBase, определённые в модуле MatrixFunctions

Остальная часть страницы документирует следующие методы MatrixBase, определённые в модуле MatrixFunctions.

MatrixBase::cos()

Вычисление матричного косинуса.

const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
Параметры
[in] M квадратная матрица.
Возвращает
выражение, представляющее \( \cos(M) \).

Эта функция вычисляет матричный косинус. Используйте ArrayBase::cos() для вычисления поэлементного косинуса.

Реализация вызывает matrixFunction() со StdStemFunctions::cos().

См. также
sin() для примера.

MatrixBase::cosh()

Вычисление матричного гиперболического косинуса.

const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const
Параметры
[in] M квадратная матрица.
Возвращает
выражение, представляющее \( \cosh(M) \)

Эта функция вызывает matrixFunction() со StdStemFunctions::cosh().

См. также
sinh() для примера.

MatrixBase::exp()

Вычисление матричной экспоненты.

const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
Параметры
[in] M матрица, экспонента которой должна быть вычислена.
Возвращает
выражение, представляющее матричную экспоненту M.

Матричная экспонента \( M \) определяется по формуле

\[ \exp(M) = \sum_{k=0}^\infty \frac{M^k}{k!}. \]

Матричная экспонента может быть использована для решения линейных обыкновенных дифференциальных уравнений: решение уравнения \( y' = My \) с начальным условием \( y(0) = y_0 \) задаётся формулой \( y(t) = \exp(M) y_0 \).

Матричная экспонента отличается от применения функции exp ко всем элементам матрицы. Используйте ArrayBase::exp(), если хотите сделать последнее.

Стоимость вычисления приблизительно \( 20 n^3 \) для матриц размера \( n \). Число 20 слабо зависит от нормы матрицы.

Матричная экспонента вычисляется с помощью метода масштабирования и возведения в квадрат в сочетании с приближением Паде. Матрица сначала масштабируется, затем вычисляется приближение экспоненты уменьшенной матрицы, а затем масштабирование отменяется с помощью повторного возведения в квадрат. Степень приближения Паде выбирается таким образом, чтобы погрешность приближения была меньше погрешности округления. Однако ошибки могут накапливаться на стадии возведения в квадрат.

Подробности алгоритма можно найти в: Nicholas J. Higham, "The scaling and squaring method for the matrix exponential revisited," SIAM J. Matrix Anal. Applic., 26:1179–1193, 2005.

Пример: следующая программа проверяет, что

\[ \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]. \]

Это соответствует вращению на \( \frac14\pi \) радиан вокруг оси z.

#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";
}

Вывод:

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 должна быть матрицей типа float, double, long double complex<float>, complex<double>, или complex<long double>.

MatrixBase::log()

Вычисление матричного логарифма.

const MatrixLogarithmReturnValue<Derived> MatrixBase<Derived>::log() const
Параметры
[in] M обратимая матрица, логарифм которой необходимо вычислить.
Возвращает
выражение, представляющее матричный логарифм корня из M.

Матричный логарифм \( M \) — это матрица \( X \), такая что \( \exp(X) = M \), где exp обозначает матричную экспоненту. Как и для скалярного логарифма, уравнение \( \exp(X) = M \) может иметь несколько решений; эта функция возвращает матрицу, у которой собственные значения имеют мнимую часть в интервале \( (-\pi,\pi] \).

Матричный логарифм отличается от применения функции log ко всем элементам матрицы. Используйте ArrayBase::log(), если хотите сделать последнее.

В вещественном случае матрица \( M \) должна быть обратимой и не должна иметь собственных значений, которые являются вещественными и отрицательными (пары комплексно-сопряжённых собственных значений разрешены). В комплексном случае она должна быть только обратимой.

Эта функция вычисляет матричный логарифм с помощью алгоритма Шура-Парлетта, реализованного в MatrixBase::matrixFunction(). Логарифм атомного блока вычисляется с помощью MatrixLogarithmAtomic, который использует прямое вычисление для блоков 1×1 и 2×2 и алгоритм обратного масштабирования и возведения в квадрат для больших блоков, при этом квадратные корни вычисляются с помощью MatrixBase::sqrt().

Подробности алгоритма можно найти в разделе 11.6.2: Nicholas J. Higham, Functions of Matrices: Theory and Computation, SIAM 2008. ISBN 978-0-898716-46-7.

Пример: следующая программа проверяет, что

\[ \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]. \]

Это соответствует вращению на \( \frac14\pi \) радиан вокруг оси z. Это обратное к примеру, используемому в документации exp().

#include <unsupported/Eigen/MatrixFunctions>
#include <iostream>
 
using namespace Eigen;
 
int main()
{
  using std::sqrt;
  MatrixXd A(3,3);
  A << 0.5*sqrt(2), -0.5*sqrt(2), 0,
       0.5*sqrt(2),  0.5*sqrt(2), 0,
       0,            0,           1;
  std::cout << "The matrix A is:\n" << A << "\n\n";
  std::cout << "The matrix logarithm of A is:\n" << A.log() << "\n";
}

Вывод:

The matrix A is:
 0.707107 -0.707107         0
 0.707107  0.707107         0
        0         0         1

The matrix logarithm of A is:
-8.86512e-17    -0.785398            0
    0.785398 -8.86512e-17            0
           0            0            0
Примечание
M должна быть матрицей типа float, double, long double, complex<float>, complex<double>, или complex<long double>.
См. также
MatrixBase::exp(), MatrixBase::matrixFunction(), класс MatrixLogarithmAtomic, MatrixBase::sqrt().

MatrixBase::pow()

Вычисление матрицы в степени произвольной действительной степени.

const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(RealScalar p) const
Параметры
[in] M основание матричной степени, должна быть квадратной матрицей.
[in] p показатель степени матричной степени.

Матричная степень \( M^p \) определяется как \( \exp(p \log(M)) \), где exp обозначает матричную экспоненту, а log обозначает матричный логарифм. Это отличается от возведения всех элементов матрицы в p-ю степень. Используйте ArrayBase::pow(), если хотите сделать последнее.

Если p является комплексным числом, тип скаляра M должен соответствовать типу p . \( M^p \) просто вычисляется как \( \exp(p \log(M)) \). Поэтому матрица \( M \) должна удовлетворять условиям, чтобы быть аргументом матричного логарифма.

Если p является вещественным числом, оно приводится к типу вещественного скаляра M. Затем эта функция вычисляет матричную степень с помощью алгоритма Шура-Паде, реализованного в классе MatrixPower. Показатель степени разбивается на целую и дробную часть, где дробная часть находится в интервале \( (-1, 1) \). Главный диагональный элемент и первый наддиагональный элемент вычисляются непосредственно.

Если M является сингулярной матрицей с полупростым собственным нулевым значением, а p является положительным, то фактор Шура \( T \) переупорядочивается с помощью вращений Гивенса, т.е.

\[ T = \left[ \begin{array}{cc} T_1 & T_2 \\ 0 & 0 \end{array} \right] \]

где \( T_1 \) обратима. Тогда \( T^p \) задаётся как

\[ T^p = \left[ \begin{array}{cc} T_1^p & T_1^{-1} T_1^p T_2 \\ 0 & 0 \end{array}. \right] \]

Предупреждение
Дробная степень матрицы с неполупростым нулевым собственным значением не определена однозначно. Мы вводим отказ от проверки против неточного результата, например,
#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;
}
.

Подробности алгоритма можно найти в: Nicholas J. Higham and Lijing Lin, "A Schur-Pad&eacute; algorithm for fractional powers of a matrix," SIAM J. Matrix Anal. Applic., 32(3):1056–1078, 2011.

Пример: следующая программа проверяет, что

\[ \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]. \]

Это соответствует вращению на \( \frac14\pi \) радиан вокруг оси z.

#include <unsupported/Eigen/MatrixFunctions>
#include <iostream>
 
using namespace Eigen;
 
int main()
{
  const double pi = std::acos(-1.0);
  Matrix3d A;
  A << cos(1), -sin(1), 0,
       sin(1),  cos(1), 0,
           0 ,      0 , 1;
  std::cout << "The matrix A is:\n" << A << "\n\n"
               "The matrix power A^(pi/4) is:\n" << A.pow(pi/4) << std::endl;
  return 0;
}

Вывод:

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() удобен для пользователя. Однако в некоторых случаях следует использовать класс MatrixPower напрямую. MatrixPower может сохранять результат разложения Шура, поэтому он лучше подходит для вычисления различных степеней для одной и той же матрицы.

Пример:

#include <unsupported/Eigen/MatrixFunctions>
#include <iostream>
 
using namespace Eigen;
 
int main()
{
  Matrix4cd A = Matrix4cd::Random();
  MatrixPower<Matrix4cd> Apow(A);
 
  std::cout << "The matrix A is:\n" << A << "\n\n"
               "A^3.1 is:\n" << Apow(3.1) << "\n\n"
               "A^3.3 is:\n" << Apow(3.3) << "\n\n"
               "A^3.7 is:\n" << Apow(3.7) << "\n\n"
               "A^3.9 is:\n" << Apow(3.9) << std::endl;
  return 0;
}

Вывод:

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 должна быть матрицей типа float, double, long double, complex<float>, complex<double>, или complex<long double>.
См. также
MatrixBase::exp(), MatrixBase::log(), класс MatrixPower.

MatrixBase::matrixFunction()

Вычисление матричной функции.

const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename internal::stem_function<typename internal::traits<Derived>::Scalar>::type f) const
Параметры
[in] M аргумент матричной функции, должна быть квадратной матрицей.
[in] f целая функция; f(x,n) должна вычислять n-ю производную f в точке x.
Возвращает
выражение, представляющее f, применённое к M.

Предположим, что M — это матрица, элементы которой имеют тип Scalar. Тогда второй аргумент, f, должен быть функцией с прототипом

ComplexScalar f(ComplexScalar, int) 

где ComplexScalar = std::complex<Scalar> если Scalar вещественное (например, float или double) и ComplexScalar = Scalar если Scalar комплексное. Значение f(x,n) должно быть \( f^{(n)}(x) \), n-й производной f в точке x.

Эта процедура использует алгоритм, описанный в: Philip Davies and Nicholas J. Higham, "A Schur-Parlett algorithm for computing matrix functions", SIAM J. Matrix Anal. Applic., 25:464–485, 2003.

Фактическая работа выполняется классом MatrixFunction.

Пример: следующая программа проверяет, что

\[ \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]. \]

Это соответствует повороту на \( \frac14\pi \) радиана вокруг оси z. Это тот же пример, что используется в документации 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";
}

Вывод:

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

Обратите внимание, что функция expfn определена для комплексных чисел x, даже если матрица A определена над вещественными числами. Вместо expfn, мы также могли использовать StdStemFunctions::exp:

A.matrixFunction(StdStemFunctions<std::complex<double> >::exp, &B);

MatrixBase::sin()

Вычислить синус матрицы.

const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
Параметры
[in] M квадратная матрица.
Возвращает
выражение, представляющее \( \sin(M) \).

Эта функция вычисляет синус матрицы. Для вычисления поэлементного синуса используйте ArrayBase::sin().

Реализация вызывает matrixFunction() со значением StdStemFunctions::sin().

Пример:

#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";
}

Вывод:

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.77156e-16  4.71845e-16
-5.55112e-17            1  2.77556e-16
 1.66533e-16 -2.08167e-16            1

MatrixBase::sinh()

Вычислить гиперболический синус матрицы.

MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const
Параметры
[in] M квадратная матрица.
Возвращает
выражение, представляющее \( \sinh(M) \)

Эта функция вызывает matrixFunction() со значением StdStemFunctions::sinh().

Пример:

#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";
}

Вывод:

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.19209e-07           0
2.10479e-07           1 1.78814e-07
1.19209e-07 2.38419e-07           1

MatrixBase::sqrt()

Вычислить квадратный корень матрицы.

const MatrixSquareRootReturnValue<Derived> MatrixBase<Derived>::sqrt() const
Параметры
[in] M обратимая матрица, для которой необходимо вычислить квадратный корень.
Возвращает
выражение, представляющее квадратный корень матрицы M.

Квадратный корень матрицы \( M \) — это матрица \( M^{1/2} \), квадрат которой равен исходной матрице; таким образом, если \( S = M^{1/2} \), то \( S^2 = M \). Это отличается от извлечения квадратного корня из всех элементов матрицы; используйте ArrayBase::sqrt(), если нужно сделать последнее.

В вещественном случае матрица \( M \) должна быть обратимой, и у неё не должно быть вещественных отрицательных собственных значений (допускаются пары комплексно сопряжённых собственных значений). В этом случае матрица имеет вещественный квадратный корень, и это тот квадратный корень, который вычисляет эта функция.

Квадратный корень матрицы вычисляется путём приведения матрицы к квазитреугольной форме с помощью вещественного разложения Шура. Квадратный корень квазитреугольной матрицы затем можно вычислить напрямую. Стоимость составляет приблизительно \( 25 n^3 \) действительных операций с плавающей точкой для вещественного разложения Шура и \( 3\frac13 n^3 \) действительных операций с плавающей точкой для остальной части (хотя время вычисления на практике, вероятно, больше, чем это показывает).

Подробности алгоритма можно найти в: Nicholas J. Highan, "Computing real square roots of a real matrix", Linear Algebra Appl., 88/89:405–430, 1987.

Если матрица положительно определённая симметричная, то квадратный корень также будет положительно определённой симметричной матрицей. В этом случае лучше использовать SelfAdjointEigenSolver::operatorSqrt() для его вычисления.

В комплексном случае матрица \( M \) должна быть обратимой; это ограничение алгоритма. Квадратный корень, вычисленный этим алгоритмом, — это тот, у которого аргументы собственных значений лежат в интервале \( (-\frac12\pi, \frac12\pi] \). Это стандартный разрез ветви.

Вычисления аналогичны вещественному случаю, за исключением того, что используется комплексное разложение Шура для приведения матрицы к треугольной форме. Теоретическая стоимость такая же. Подробности см.: Åke Björck и Sven Hammarling, "A Schur method for the square root of a matrix", Linear Algebra Appl., 52/53:127–140, 1983.

Пример: следующая программа проверяет, что квадратный корень

\[ \left[ \begin{array}{cc} \cos(\frac13\pi) & -\sin(\frac13\pi) \\ \sin(\frac13\pi) & \cos(\frac13\pi) \end{array} \right], \]

соответствующий повороту на 60 градусов, является поворотом на 30 градусов:

\[ \left[ \begin{array}{cc} \cos(\frac16\pi) & -\sin(\frac16\pi) \\ \sin(\frac16\pi) & \cos(\frac16\pi) \end{array} \right]. \]

#include <unsupported/Eigen/MatrixFunctions>
#include <iostream>
 
using namespace Eigen;
 
int main()
{
  const double pi = std::acos(-1.0);
 
  MatrixXd A(2,2);
  A << cos(pi/3), -sin(pi/3), 
       sin(pi/3),  cos(pi/3);
  std::cout << "The matrix A is:\n" << A << "\n\n";
  std::cout << "The matrix square root of A is:\n" << A.sqrt() << "\n\n";
  std::cout << "The square of the last matrix is:\n" << A.sqrt() * A.sqrt() << "\n";
}

Вывод:

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
См. также
класс RealSchur, класс ComplexSchur, класс MatrixSquareRoot, SelfAdjointEigenSolver::operatorSqrt().

matrix_sqrt_quasi_triangular()

template<typename MatrixType , typename ResultType >
void Eigen::matrix_sqrt_quasi_triangular ( const MatrixType & arg,
ResultType & result
)

Вычислить квадратный корень квазитреугольной матрицы.

Шаблонные параметры
MatrixType тип arg, аргумента квадратного корня матрицы, ожидается, что это будет экземпляр шаблона класса Matrix.
ResultType тип result, в котором будет храниться результат.
Параметры
[in] arg аргумент квадратного корня матрицы.
[out] result квадратный корень верхней квазитреугольной части arg.

Эта функция вычисляет квадратный корень верхней квазитреугольной матрицы, хранящейся в верхней квазитреугольной части arg. Только верхняя квазитреугольная часть result обновляется, остальная часть не затрагивается. Подробности о том, как реализовано это вычисление, см. в MatrixBase::sqrt().

См. также
MatrixSquareRoot, MatrixSquareRootQuasiTriangular

matrix_sqrt_triangular()

template<typename MatrixType , typename ResultType >
void Eigen::matrix_sqrt_triangular ( const MatrixType & arg,
ResultType & result
)

Вычислить квадратный корень треугольной матрицы.

Шаблонные параметры
MatrixType тип arg, аргумента квадратного корня матрицы, ожидается, что это будет экземпляр шаблона класса Matrix.
ResultType тип result, в котором будет храниться результат.
Параметры
[in] arg аргумент квадратного корня матрицы.
[out] result квадратный корень верхней треугольной части arg.

Обновляется только верхняя треугольная часть (включая диагональ) result, остальная часть не затрагивается. Подробности о том, как реализовано это вычисление, см. в MatrixBase::sqrt().

См. также
MatrixSquareRoot, MatrixSquareRootQuasiTriangular

© Eigen.
Licensed under the MPL2 License.
https://eigen.tuxfamily.org/dox/unsupported/group__MatrixFunctions__Module.html

Spec-Zone.ru

Настройки Оффлайн Что нового Помощь О нас
Spec-Zone .ru
спецификации, руководства, описания, API