С ++ Холецкая факторизация

Question

Мне нужно переписать код MatLab с использованием C ++.

Внутри кода Matlab мы вызываем функцию chol для вычисления верхней матрицы tri angular.

Для части C ++ я начинаю смотреть на Eigen. Однако я изо всех сил пытаюсь получить эквивалент функции Matlab chol.

Я пытался использовать класс Eigen LDLT:

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;

int main() {

  MatrixXd matA(2, 2);
  matA << 1, 2, 3, 4;

  MatrixXd matB(4, 4);
  matB << matA, matA/10, matA/10, matA;
  matB = matB*matB.transpose();

  Eigen::LDLT<MatrixXd> tmp(matB);
  MatrixXd U = tmp.matrixU();
  cout << U << endl;
}

, но результат отличается от кода Matlab:

matB = [  1   2 0.1 0.2
          3   4 0.3 0.4
        0.1 0.2   1   2
        0.3 0.4   3   4];
matB = matB*matB';
D = chol(matB);

Answer 1

Используя ваш пример кода и документацию Matlab , я получаю тот же результат при использовании LLT вместо LDLT ( онлайн ):

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using std::cout;

int main()
{
  MatrixXd matA(3,3);
  matA << 1, 0, 1, 0, 2, 0, 1, 0, 3;
  cout << matA << "\n\n";
  Eigen::LDLT<MatrixXd> tmp(matA);
  cout << ((tmp.info() == Success) ? "succeeded" : "failed") << "\n\n";
  MatrixXd U = tmp.matrixL();
  cout << U << "\n\n";
  // Using LLT instead
  cout << MatrixXd(matA.llt().matrixL()) << "\n\n";
  cout << MatrixXd(matA.llt().matrixU()) << "\n\n";
}

Выходы:

1 0 1
0 2 0
1 0 3

успешно

1 0 0
0 1 0
0,333333 0 1

1 0 0
0 1,41421 0
1 0 1,41421

1 0 1
0 1,41421 0
0 0 1,41421

Answer 2

Согласно документации Eigen :: LDTL 2-й параметр шаблона _UpLo по умолчанию равен Lower, но вы пропустили этот параметр и хотите вычислить верхнюю матрицу tri angular.

Таким образом, ваш экземпляр класса должен выглядеть примерно так (не знаю, правильное ли здесь определение Eigen Upper):

Eigen::LDLT<MatrixXd, Upper> tmp(matB);

Edit @chtz верно - использование Upper не даст ожидаемого результата, поскольку класс LDLT предназначен для надежного разложения по Холески с поворотом. Таким образом, в дополнение к правильному ответу @Avi вы также можете использовать правильный класс для стандартного разложения по Холецкому:

Eigen::LLT<MatrixXd> tmp(matA);
cout <<  MatrixXd(tmp.matrixU()) << "\n\n";