LU разложение с PartialPivLU

У меня есть следующий код в MatLab, который выполняет декомпозицию LU, прежде чем пытаться вычислить «барашек», который я включил, чтобы дать некоторый контекст.

P=[1,2,3;4,5,6;7,8,9];
U=[0;1;2];
[F,J]=lu(P);
Jlamda=F\U;
lamb=J\Jlamda;

F это:

0,142857142857143 1 0
0,571428571428571 0,500000000000000 1
1 0 0

U это:

7 8 9
0 0.857142857142857 1.71428571428571
0 0 1.11022302462516e-16

Когда я пытаюсь повторить это в Eigen с помощью следующего кода:

MatrixXd P(3, 3);
P << 1, 2, 3, 4, 5, 6, 7, 8, 9;
MatrixXd U(3, 1);
U << 0, 1, 2;

PartialPivLU<MatrixXd> lu = PartialPivLU<MatrixXd>(P);
MatrixXd J = lu.matrixLU().triangularView<UpLoType::Upper>();
MatrixXd F = lu.matrixLU().triangularView<UpLoType::UnitLower>();

MatrixXd Jlamda = F.lu().solve(U);
MatrixXd l = J.lu().solve(Jlamda);

cout << F << endl;
cout << endl;
cout << J << endl;

Какие отпечатки:

1 0 0
0,142857 1 0
0,571429 0,5 1
7 8 9
0 0,857143 1,71429
0 0 1.11022e-16

Хотя я, очевидно, могу изготовить матрицу для преобразования строк F в C ++ в матрицу F из MatLab, я не уверен, как это сделать динамически.

Является ли PartialPivLU лучшим способом для этого, или я упускаю что-то более тривиальное?

Заранее спасибо.

1

Решение

По телефону [F,J]=lu(P)полученная матрица F это переставляются нижняя треугольная матрица. Вы можете вызвать функцию как [F,J,perm]=lu(P) получить F как действительно более низкая матрица треугольника и P в качестве отдельной матрицы перестановок, так что F*J = perm*P, По умолчанию строка

MatrixXd F = lu.matrixLU().triangularView<UpLoType::UnitLower>();

использование Eigen возвращает истинную нижнюю треугольную матрицу. Если вы хотите переставить нижнюю треугольную матрицу, такую ​​как Matlab, вы можете сохранить матрицу перестановок в Eigen, вызвав permutationP а затем умножить эту матрицу на F,

4

Другие решения

Ответ DCSmith довольно близок, однако, кажется, что вы также должны вызвать transpose () для матрицы перестановок, чтобы получить правильный результат:

MatrixXd P(3, 3);
P << 1, 2, 3, 4, 5, 6, 7, 8, 9;
MatrixXd U(3, 1);
U << 0, 1, 2;

PartialPivLU<MatrixXd> lu = PartialPivLU<MatrixXd>(P);
MatrixXd J = lu.matrixLU().triangularView<UpLoType::Upper>();
MatrixXd F = lu.matrixLU().triangularView<UpLoType::UnitLower>();
MatrixXd perm = lu.permutationP().transpose();
MatrixXd F1 = perm * F;

MatrixXd Jlamda = F.lu().solve(U);
MatrixXd l = J.lu().solve(Jlamda);

cout << perm << endl;
cout << endl;
cout << F1 << endl;

Печать:

0 1 0
0 0 1
1 0 0
0,142857 1 0
0,571429 0,5 1
1 0 0

Что совпадает с примером MatLab

1

По вопросам рекламы [email protected]