如何将特征库用于LU分解C

how to use eigen library for lu decomposition c++

本文关键字:LU 分解 用于 特征      更新时间:2023-10-16

我正在使用标准特征库来计算给定矩阵的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();
}

分解是 PA=LU,这意味着产品 LU代表矩阵A,在某些行互换之后。这与您所有三个示例一致。有关更多详细信息,请参见各自的文档。