如何使用特徵庫進行lu分解C++

Question

我正在使用標準特徵庫來計算給定矩陣的LU分解。

但我對某些功能感到困惑。問題是基於將A更改爲LU（A = LU）。

我以爲我知道LU分解背後的邏輯。首先檢查你是否可以直接做，如果不是原來的矩陣A改爲PA並做PA=LU;我認爲這一步是通過調用函數matrixLU()完成的，這正是實例2中所做的。 但是，這個函數給了我一些我在示例1和3中不能理解的東西。

示例1的作品，問題是函數matrixLU（）的用法是什麼？

例2的作品，我知道在這裏使用matrixLU()。 （更改爲PA = LU）

示例3不起作用，結果LU給出了一個矩陣，它改變了原始矩陣的行。爲什麼這個產生錯誤的結果？

欣賞函數matrixLU（）背後的任何解釋或數學邏輯。

#include <iostream> 
#include <Eigen/Dense> 
//#include <Dense> 
#include <Eigen/LU> 
using Eigen::MatrixXd; 
using std::cout; 
using std::endl; 
using Eigen::VectorXd; 
using std::cin; 
int main() 
{ 

    //example 1 
    typedef Eigen::Matrix<double,4,4> M4x4; 
    M4x4 p; 
    p<< 7,3,-1,2,3,8,1,-1,-1,1,4,-1,2,-4,-1,6; 
    cout<<p<<endl<<endl; 
    // Create LU Decomposition template object for p 
    Eigen::PartialPivLU<M4x4> LU(p); 
    cout<<"LU MATRIX:\n"<<LU.matrixLU()<<endl<<endl; 
    // Output L, the lower triangular matrix 
    M4x4 l = MatrixXd::Identity(4,4);//默認 單位對角矩陣 
    //開始填充 
    l.block<4,4>(0,0).triangularView<Eigen::StrictlyLower>()=LU.matrixLU(); 
    cout<<"L MATRIX:\n"<<l<<endl <<endl; 
    M4x4 u = LU.matrixLU().triangularView<Eigen::Upper>(); 
    cout<<"R MATRIX:\n"<<u<<endl<<endl; 
    MatrixXd m0(4,4); 
    m0 = l*u; 
    cout<<"calculate the original matrix:\n"<<m0<<endl<<endl; 
    //證明完畢 
    //example 2 
    typedef Eigen::Matrix<double,2,2> M2X2; 
    M2X2 p0; 
    p0<< 0,2,1,3; 
    cout<<p0<<endl<<endl; 
    Eigen::PartialPivLU<M2X2> LU0(p0); 
    cout<<"LU MATRIX:\n"<<LU0.matrixLU()<<endl<<endl;//原來是在做PA的過程 
    //一切結果從PA開始 
    M2X2 l0 = MatrixXd::Identity(2,2); 
    l0.block<2,2>(0,0).triangularView<Eigen::StrictlyLower>()=LU0.matrixLU(); 
    cout<<"L MATRIX:\n"<<l0<<endl <<endl; 
    //以下省略N行 
    //example 3 
    typedef Eigen::Matrix<double,3,3> M3X3; 
    M3X3 p1; 
    p1<<3,-1,2,6,-1,5,-9,7,3; 
    cout<<p1<<endl<<endl; 
    Eigen::PartialPivLU<M3X3> LU1(p1); 
    cout<<"LU MATRIX:\n"<<LU1.matrixLU()<<endl<<endl;//暫時沒明白這步做的啥 
    M3X3 l1 = MatrixXd::Identity(3,3); 
    l1.block<3,3>(0,0).triangularView<Eigen::StrictlyLower>()=LU1.matrixLU(); 
    cout<<"L MATRIX:\n"<<l1<<endl <<endl; 
    //直接up 
    M3X3 u1 = LU1.matrixLU().triangularView<Eigen::Upper>(); 
    cout<<"R MATRIX:\n"<<u1<<endl<<endl; 
    cout<<l1*u1<<endl; 
    cin.get(); 
} 

Answer 1

的因式分解是PA=LU指該產品LU表示一些行交換後的矩陣A。這與你所有的三個例子一致。有關更多詳細信息，請參閱相應的documentation。