Eigen不同的方法来求矩阵的逆的效率

Feisy

已于 2022-05-27 15:41:51 修改

阅读量4.8k

点赞数 1

分类专栏： C++ 文章标签：矩阵线性代数

于 2022-05-27 14:42:21 首次发布

本文链接：https://blog.csdn.net/feisy/article/details/125003318

版权

C++ 专栏收录该内容

24 篇文章 0 订阅

订阅专栏

背景

不同尺寸的矩阵，求逆使用不同的方法，会有不同的效率的.

16x16矩阵的直接求逆与PartialPivLU的效率对比(本人亲测）

1 实现代码

GetSystemTimeInMacroSecond的实现

    boost::posix_time::ptime now = boost::posix_time::microsec_clock::universal_time() + boost::posix_time::hours(8);
    long long cur_time = now.time_of_day().total_microseconds();
    return cur_time;

求逆用时评估代码

Eigen::MatrixXf x;
    Eigen::MatrixXf A = Eigen::MatrixXf::Zero(22, 22);
    Eigen::VectorXf b = Eigen::MatrixXf::Zero(22, 1);
    A << 1.34600000000000, 1.13841228110978e-18, -3.54804631550289e-17, -2.05829006182795e-19, -0.160432023690557, 0.00411124458769287, -0.000223899921724580, -0.0593112659337829, -0.0241839767166998, -0.00834748033971202, -5.37814422485606e-06, 5.23214930710842e-06, -0.0580682816915559, -0.0247053374485193, -0.00834723095360454, 1.97819241321089e-07, 0.999999990195040, 9.97072788916781e-05, 9.83279179330304e-05, 0, 0, 0,
    1.13841228110978e-18, 1.34600000000000, -3.32036915323686e-18, 0.160484535740467, 5.65533818549124e-05, 0.0112007248696078, 0.0837289215044391, -9.20550154007569e-05, 4.95864766071101e-05, -2.26259635075635e-05, 0.00198399199113416, 0.0818535878168905, 0.000219747293578020, -9.46046009141372e-05, 8.32360721203644e-07, 0.00198399998046782, -9.97169861081831e-05, 0.999999990155153, 9.87229323286820e-05, 0, 0, 0,
    - 3.54804631550289e-17, -3.32036915323686e-18, 1.34600000000000, 0.000324983605022713, -0.00204205102000001, 7.50143433983415e-05, -0.00174330315904512, 0.00320107672037773, 0.00538664534138679, -2.49931385726132e-06, 1.68957988594421e-06, 0.00221636220236278, -0.00770414705602320, 0.00350024291275525, 8.20683675049907e-07, -1.95885748867475e-07, -9.83180735701035e-05, -9.87327363243321e-05, 0.999999990292701, 0, 0, 0,
    - 2.05829006182795e-19, 0.160484535740467, 0.000324983605022713, 0.0368009077713332, 0.000450394440556140, 0.00224060694237125, 0.0171499333910333, 0.000111460617639351, 0.000231822191509620, -6.19505197893729e-06, 0.000604514119035474, 0.0167286230613098, 0.000356657959376959, -0.000161229251681464, 1.92165660849536e-07, 0.000604838471732968, -2.34113026654352e-05, 0.273495197753565, -0.0392376814557259, 0.999999990195040, 9.97072788916781e-05, 9.83279179330304e-05,
    - 0.160432023690557, 5.65533818549124e-05, -0.00204205102000001, 0.000450394440556140, 0.0349927346271823, -0.000770239409681646, -5.52580627558622e-06, 0.0135113477743587, 0.00619537343699084, 0.00225656479992005, 6.33031697594492e-07, -2.32778930395674e-07, 0.0132115377315313, 0.00612436104238307, 0.00223967166698827, 5.63540677387556e-07, -0.272298424113141, 0.000287254259750317, -0.0113390517979790, -0.000102400584293418, 0.999613987367637, 0.0277824724078660,
    0.00411124458769287, 0.0112007248696078, 7.50143433983415e-05, 0.00224060694237125, -0.000770239409681646, 0.00250513838063355, -7.96998049017821e-05, 0.00207932171781462, 0.000827490251987481, 0.000276737538215227, -2.91526011283910e-06, 0.00218079218578034, -0.00265933192175922, -0.00115248022548862, -0.000388572073708229, 5.69402771431904e-05, 0.0467420522622825, 0.0269122305302098, -0.00192157718616087, 0.0569942750587246, -0.0277314782938061, 0.997989287378462,
    - 0.000223899921724580, 0.0837289215044391, -0.00174330315904512, 0.0171499333910333, -5.52580627558622e-06, -7.96998049017821e-05, 0.0173695770953092, 9.19953701202703e-08, 5.55632170658309e-07, -2.45754961919456e-07, 0.000608127257267949, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    - 0.0593112659337829, -9.20550154007569e-05, 0.00320107672037773, 0.000111460617639351, 0.0135113477743587, 0.00207932171781462, 9.19953701202703e-08, 0.0107948324950528, 0.00507773016911412, 0.00187572295300195, 7.24739201536356e-07, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    - 0.0241839767166998, 4.95864766071101e-05, 0.00538664534138679, 0.000231822191509620, 0.00619537343699084, 0.000827490251987481, 5.55632170658309e-07, 0.00507773016911412, 0.00297121709317547, 0.00113838089040340, 9.65602255182848e-07, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    - 0.00834748033971202, -2.26259635075635e-05, -2.49931385726132e-06, -6.19505197893729e-06, 0.00225656479992005, 0.000276737538215227, -2.45754961919456e-07, 0.00187572295300195, 0.00113838089040340, 0.000601080937631330, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    - 5.37814422485606e-06, 0.00198399199113416, 1.68957988594421e-06, 0.000604514119035474, 6.33031697594492e-07, -2.91526011283910e-06, 0.000608127257267949, 7.24739201536356e-07, 9.65602255182848e-07, 0, 0.000121744000000000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    5.23214930710842e-06, 0.0818535878168905, 0.00221636220236278, 0.0167286230613098, -2.32778930395674e-07, 0.00218079218578034, 0, 0, 0, 0, 0, 0.0166289413041880, 1.54987609437772e-06, -2.68840504135592e-07, 3.06177337509682e-07, 0.000600039585916120, -4.76000849622613e-05, 0.271127384428256, 0.00671673095904333, 0.998374505182241, 0.00158666949874198, -0.0569721851609907,
    - 0.0580682816915559, 0.000219747293578020, -0.00770414705602320, 0.000356657959376959, 0.0132115377315313, -0.00265933192175922, 0, 0, 0, 0, 0, 1.54987609437772e-06, 0.0105650430359955, 0.00498056352012310, 0.00185632506252512, 7.71547978523958e-07, -0.226373410990567, 0.000388885069820543, -0.0139636088181738, 0, 0.999612417060891, 0.0278391030330401,
    - 0.0247053374485193, -9.46046009141372e-05, 0.00350024291275525, -0.000161229251681464, 0.00612436104238307, -0.00115248022548862, 0, 0, 0, 0, 0, -2.68840504135592e-07, 0.00498056352012310, 0.00300667325425074, 0.00115608170315948, -7.85620278703043e-07, -0.142484112085114, -0.000395977963055969, 0.0142182917417796, 0, 0.999612417060891, 0.0278391030330401,
    - 0.00834723095360454, 8.32360721203644e-07, 8.20683675049907e-07, 1.92165660849536e-07, 0.00223967166698827, -0.000388572073708229, 0, 0, 0, 0, 0, 3.06177337509682e-07, 0.00185632506252512, 0.00115608170315948, 0.000601026402068215, 0, -0.0759883725118268, 0, 0, 0, 0.999612417060891, 0.0278391030330401,
    1.97819241321089e-07, 0.00198399998046782, -1.95885748867475e-07, 0.000604838471732968, 5.63540677387556e-07, 5.69402771431904e-05, 0, 0, 0, 0, 0, 0.000600039585916120, 7.71547978523958e-07, -7.85620278703043e-07, 0, 0.000121744000000000, 0, 0.0300000000000000, 0, 1, 0, 0,
    0.999999990195040, -9.97169861081831e-05, -9.83180735701035e-05, -2.34113026654352e-05, -0.272298424113141, 0.0467420522622825, 0, 0, 0, 0, 0, -4.76000849622613e-05, -0.226373410990567, -0.142484112085114, -0.0759883725118268, 0, 0, 0, 0, 0, 0, 0,
    9.97072788916781e-05, 0.999999990155153, -9.87327363243321e-05, 0.273495197753565, 0.000287254259750317, 0.0269122305302098, 0, 0, 0, 0, 0, 0.271127384428256, 0.000388885069820543, -0.000395977963055969, 0, 0.0300000000000000, 0, 0, 0, 0, 0, 0,
    9.83279179330304e-05, 9.87229323286820e-05, 0.999999990292701, -0.0392376814557259, -0.0113390517979790, -0.00192157718616087, 0, 0, 0, 0, 0, 0.00671673095904333, -0.0139636088181738, 0.0142182917417796, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0.999999990195040, -0.000102400584293418, 0.0569942750587246, 0, 0, 0, 0, 0, 0.998374505182241, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 9.97072788916781e-05, 0.999613987367637, -0.0277314782938061, 0, 0, 0, 0, 0, 0.00158666949874198, 0.999612417060891, 0.999612417060891, 0.999612417060891, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 9.83279179330304e-05, 0.0277824724078660, 0.997989287378462, 0, 0, 0, 0, 0, -0.0569721851609907, 0.0278391030330401, 0.0278391030330401, 0.0278391030330401, 0, 0, 0, 0, 0, 0, 0;


    b << 0,
    0,
    -13.2042600000000,
    -0.00318808916527274,
    0.0200325205062001,
    -0.000735890708737724,
    0.0155774925360128,
    -0.0369437528953341,
    -0.0409522547984081,
    -0.00390667008137935,
    -0.00232355145046985,
    0.413165511114724,
    2.58954781196455,
    -0.0402923712348155,
    0.00438897429704739,
    0.000278169226818484,
    -0.0678455928345908,
    0.0527586395813534,
    -0.726306425617754,
    0.000161368614991371,
    -0.0112377101793276,
    0.237950558222788;

    long long end_time ;
    long long elapsed_time;
    long long start_time = TimeUtil::GetSystemTimeInMacroSecond();
    //
//    Eigen::Vector3f x = A.colPivHouseholderQr().solve(b);
//    std::cout << "The solution is:\n" << x << std::endl;
//    long long end_time = TimeUtil::GetSystemTimeInMacroSecond();
//    long long elapsed_time = end_time - start_time;
//    cout << "inverrse matrix use colPivHouseholderQr time = " << elapsed_time << endl;

    int test_count = 1;
    long long use_time = 0;
    for (int i = 0; i < test_count; ++i)
    {
        start_time = TimeUtil::GetSystemTimeInMacroSecond();
        Eigen::PartialPivLU<Eigen::MatrixXf> ALU(A); // A的LU分解
        Eigen::MatrixXf x1 = ALU.inverse() * b; // solve Ax=b, same as x = lu.solve(b);
        //Eigen::MatrixXd x = b * lu.inverse(); // solve xA=b
        end_time = TimeUtil::GetSystemTimeInMacroSecond();
        elapsed_time = end_time - start_time;
        use_time = use_time + elapsed_time;
    }

    cout << "inverrse matrix use PartialPivLU time = " << use_time / test_count << endl;


    use_time = 0;
    for (int i = 0; i < test_count; ++i)
    {
        start_time = TimeUtil::GetSystemTimeInMacroSecond();
        Eigen::MatrixXf x2 = A.inverse() * b;
        end_time = TimeUtil::GetSystemTimeInMacroSecond();
        elapsed_time = end_time - start_time;
        use_time = use_time + elapsed_time;
    }
    cout << "inverrse matrix use inverse time = " << use_time / test_count << endl;

2 测试环境:

OS:win11
硬件：Intel® Core™ i7-10510U CPU @ 1.80GHz 2.30 GHz
开发工具qt+vs编译器

3 debug计算时间（注意，时间跟配置和操作系统很大关系，个人不同，这里只是相对比较）

这里不是确定的值是因为，我会执行多次，观察值所在范围

只执行一次：变动很大，0-1000微秒都出现过
执行100次：PartialPivLU用时500-600多微秒，直接求逆用时300-400维秒
执行1000次：PartialPivLU用时300多微秒，直接求逆用时200-300多维秒
执行10000次：PartialPivLU用时[300,350]微秒，直接求逆用时[300,350]维秒

4 release计算时间（注意，时间跟配置和操作系统很大关系，个人不同，这里只是相对比较）

这里不是确定的值是因为，我会执行多次，观察值所在范围

只执行一次：变动很大，0微秒
执行100次：PartialPivLU用时[10,40]微秒，直接求逆用时[10,30]维秒
执行1000次：PartialPivLU用时[20,30]微秒，直接求逆用时[20,30]多维秒
执行10000次：PartialPivLU用时[15,18]微秒，直接求逆用时[12,14]维秒

4x4矩阵的求逆测试讨论

https://stackoverflow.com/questions/50909385/eigen-linear-solver-for-very-small-square-matrix

矩阵的逆的问题，一般类似求解：Ax=b----------------->x = A^-1b

这里给出stackoverflow的一个关于这个求解的的效率的讨论：

1 测试的方法以及种类

I am using Eigen on a C++ program for solving linear equation for very small square matrix(4X4).
My test code is like

template<template <typename MatrixType> typename EigenSolver>
Vertor3d solve(){
   //Solve Ax = b and A is a real symmetric matrix and positive semidefinite

   ... // Construct 4X4 square matrix A and 4X1 vector b

   EigenSolver<Matrix4d> solver(A);
   auto x = solver.solve(b);

   ... // Compute relative error for validating
}

I test some EigenSolver which include:（测试的方法种类）

FullPixLU
PartialPivLU
HouseholderQR
ColPivHouseholderQR
ColPivHouseholderQR
CompleteOrthogonalDecomposition
LDLT
Direct Inverse

Direct Inverse is:(直接求逆是调用Eigen的接口Inverse来计算逆矩阵，然后乘上去）

template<typename MatrixType>
struct InverseSolve
{
private:
    MatrixType  inv;
public:
    InverseSolve(const MatrixType &matrix) :inv(matrix.inverse()) {
    }
    template<typename VectorType>
    auto solve(const VectorType & b) {
        return inv * b;
    }
};

2 测试的数据

all use 1000000 matrices with random double from uniform distribution [0,100].I fristly construct upper-triangle and then copy to lower-triangle.

3 测试的结果：直接求逆最快，但是精度稍微差一点

I found that the fast method is DirectInverse,Even If I linked Eigen with MKL , the result was not change.

This is the test result：
FullPixLU : 477 ms
PartialPivLU : 468 ms
HouseholderQR : 849 ms
ColPivHouseholderQR : 766 ms
ColPivHouseholderQR : 857 ms
CompleteOrthogonalDecomposition : 832 ms
LDLT : 477 ms
Direct Inverse : 88 ms

The only problem of Direct Inverse is that its relative error slightly larger than other solver but acceptble.
Is there any faster or more felegant solution for my program?Is DirectI nverse the fast solution for my program?
Direct Inverse does not use the symmetric infomation so why is Direct Inverse far faster than LDLT?

For fixed sized matrices up to 4x4, inverse() is computed directly using co-factors, 
which is very fast and requires only on division, but is expected to be slightly less accurate, in general.
The solve as you implemented it, will only require a very fast Matrix-Vector product.
You can increase the accuracy by doing a refinement-step: x=inv*b; x+=inv*(b-A*x);. 
MKL will generally not help for small fixed-sized inputs

回答一：微小矩阵直接求逆比用求解法更快

Despite what many people suggest of never explicitly computing an inverse when you only want to solve a linear system, for very small matrices this can actually be beneficial, since there are closed-form solutions using co-factors.

All other alternatives you tested will be slower, since they will do pivoting (which implies branching), even for small fixed-sized matrices. Also, most of them will result in more divisions and be not vectorizable as good, as the direct computation.

To increase the accuracy (this technique can actually be used independent of the solver if required), you can refine an initial solution by solving the system again with the residual:

Eigen::Vector4d solveDirect(const Eigen::Matrix4d& A, const Eigen::Vector4d& b)
{
    Eigen::Matrix4d inv = A.inverse();
    Eigen::Vector4d x = inv * b;
    x += inv*(b-A*x);
    return x;
}

超大矩阵的求逆

https://stackoverflow.com/questions/39190245/eigen-get-inverse-failed-and-inverse-is-so-slow

Full-pivoting LU is known to be very slow, regardless of its implementation.

Better use PartialPivLU, which benefits from high performance matrix-matrix operations. Then to get the best of Eigen, use the 3.3-beta2 release and compile with both FMA (-mfma) and OpenMP (e.g., -fopenmp) supports, and don’t forget to enable compiler optimizations -O3. This operation should not take more than a few seconds.

Finally, do you really need to explicitly compute the inverse? If you only apply it to some vectors or matrices (i.e., A^-1 * B or B * A^-1) then better apply the inverse in factorized form rather than explicitly computing it. With Eigen 3.3:

MatrixXd A = ...;
PartialPivLU<MatrixXd> lu(A);
x = lu.inverse() * b; // solve Ax=b, same as x = lu.solve(b);
x = b * lu.inverse(); // solve xA=b

n these expressions, the inverse is not explicitly computed!