Внутриматричные разложения матриц

Начиная с Eigen 3.3, разложения LU, Cholesky и QR могут работать внутри заданной входной матрицы, то есть напрямую внутри неё. Эта функция особенно полезна при работе с огромными матрицами и/или когда доступная память очень ограничена (встраиваемые системы).

Для этого соответствующий класс разложения должен быть инициализирован с типом матрицы Ref<>, а объект разложения должен быть создан со входной матрицей в качестве аргумента. Рассмотрим, например, внутриматричное разложение LU с частичным выбором опор.

Начнём с основных включений и объявления 2x2 матрицы A:

#include <iostream>
#include <Eigen/Dense>
 
using namespace std;
using namespace Eigen;
 
int main()
{
  MatrixXd A(2,2);
  A << 2, -1, 1, 3;
  cout << "Here is the input matrix A before decomposition:\n" << A << endl;

Here is the input matrix A before decomposition:
 2 -1
 1  3

Ничего удивительного! Затем объявим наш объект внутриматричного LU lu, и проверим содержимое матрицы A:

  PartialPivLU<Ref<MatrixXd> > lu(A);
  cout << "Here is the input matrix A after decomposition:\n" << A << endl;

Here is the input matrix A after decomposition:
  2  -1
0.5 3.5

Здесь объект lu вычисляет и хранит факторы L и U внутри памяти, занимаемой матрицей A. Коэффициенты A были уничтожены во время факторизации и заменены факторами L и U, как можно проверить:

  cout << "Here is the matrix storing the L and U factors:\n" << lu.matrixLU() << endl;

Here is the matrix storing the L and U factors:
  2  -1
0.5 3.5

Затем можно использовать объект lu как обычно, например, для решения задачи Ax=b:

  MatrixXd A0(2,2); A0 << 2, -1, 1, 3;
  VectorXd b(2);    b << 1, 2;
  VectorXd x = lu.solve(b);
  cout << "Residual: " << (A0 * x - b).norm() << endl;

Residual: 0

Здесь, так как содержимое исходной матрицы A было потеряно, нам пришлось объявить новую матрицу A0 для проверки результата.

Поскольку память разделяется между A и lu, изменение матрицы A сделает lu недействительным. Это легко проверить, изменив содержимое A и попытавшись снова решить исходную задачу:

  A << 3, 4, -2, 1;
  x = lu.solve(b);
  cout << "Residual: " << (A0 * x - b).norm() << endl;

Residual: 15.8114

Обратите внимание, что в подкапотной части нет умных указателей, поэтому ответственность пользователя заключается в том, чтобы исходная матрица A оставалась живой, пока жив объект lu.

Если требуется обновить факторизацию изменённой A, нужно вызвать метод compute как обычно:

  A0 = A; // save A
  lu.compute(A);
  x = lu.solve(b);
  cout << "Residual: " << (A0 * x - b).norm() << endl;

Residual: 0

Обратите внимание, что вызов метода compute не изменяет память, на которую ссылается объект lu. Таким образом, если метод compute вызывается с другой матрицей A1, отличной от A, содержимое A1 не изменится. Это по-прежнему содержимое A, которое будет использовано для хранения факторов L и U матрицы A1. Это легко проверить следующим образом:

  MatrixXd A1(2,2);
  A1 << 5,-2,3,4;
  lu.compute(A1);
  cout << "Here is the input matrix A1 after decomposition:\n" << A1 << endl;

Here is the input matrix A1 after decomposition:
 5 -2
 3  4

Матрица A1 не изменена, и можно решить A1*x=b и напрямую проверить остаток без копирования A1:

  x = lu.solve(b);
  cout << "Residual: " << (A1 * x - b).norm() << endl;

Residual: 2.48253e-16

Вот список матричных разложений, поддерживающих этот механизм внутриматричных вычислений:

• класс LLT
• класс LDLT
• класс PartialPivLU
• класс FullPivLU
• класс HouseholderQR
• класс ColPivHouseholderQR
• класс FullPivHouseholderQR
• класс CompleteOrthogonalDecomposition