线性方程理论说明和Eigen解线性方程求解方法汇总

Posted 2022-02-22 非晚非晚

文章目录

• 1. 线性方程组(矩阵方程)理论
• 2. QR分解
• • 2.1 HouseholderQR
  • 2.2 ColPivHouseholderQR
  • 2.3 FullPivHouseholderQR
• 3. LLT分解
• 4. LDLT分解
• 5. LU分解
• • 5.1 LU分解
  • 5.2 partialPivLu
  • 5.3 fullPivLu
• 6. SVD分解
• • 6.1 BDCSVD分解
  • 6.2 JacobiSVD分解
• 7. 逆矩阵法

1. 线性方程组(矩阵方程)理论

线性方程组是各个方程关于未知量均为一次的方程组，方程表示如下：

$$\{\begin{array}{c}{a}_1{x}_1+b_1{x}_2+c_1{x}_3=d_1\\ a_2{x}_1+b_2{x}_2+c_2{x}_3=d_2\\ a_3{x}_1+b_3{x}_2+c_3{x}_3=d_3\end{array}$$
把线程方程写成矩阵形式如下：
$$Ax=b$$
其中，
$$A=\left(\begin{array}{ccc}{a}_1& b_1& c_1\\ a_2& b_2& c_2\\ a_3& b_3& c_3\end{array}\right),x=\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right),d=\left(\begin{array}{c}{d}_1\\ d_2\\ d_3\end{array}\right)$$

Eigen中提供了丰富的线性方程求解方法，包括LU分解法，QR分解法，SVD（奇异值分解）、特征值分解等根据A的矩阵类型、对结果的精度要求以及计算速度的要求，可以选择不同的计算方式，下面一一进行介绍。

下图是Eigen一些分解方法的简介。

如果对矩阵很了解，那么可以很方便的选择一种合理的求解方法，比如如果矩阵A是满秩且非对称矩阵，那么可以选用PartialPivLU求解方法，如果知道你的矩阵是对称正定的，那么选用LLT或者LDLT分解是一个很好的选择。下图是一些求解方法的速度。

2. QR分解

QR（正交三角）分解法是求一般矩阵全部特征值的最有效并广泛应用的方法，一般矩阵先经过正交相似变化成为Hessenberg矩阵，然后再应用QR方法求特征值和特征向量。它是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R，所以称为QR分解法，与此正规正交矩阵的通用符号Q有关。

$$A=QR$$
其中Q是正交矩阵（或酉矩阵），即$Q{Q}^T=1$，$R$是上三角矩阵。QR分解有三种常用方法：Givens 变换、Householder 变换，以及 Gram-Schmidt正交化。

• HouseholderQR：无旋转（no pivoting），速度很快但不稳定
• ColPivHouseholderQR：列旋转（column pivoting），速度稍慢但更精确
• FullPivHouseholderQR：全旋转（full pivoting），速度慢，最稳定

下面的代码无法比较速度和精度，因为矩阵太小了。

2.1 HouseholderQR

HouseholderQR分解是3种QR分解中速度最快的一种，但精度最低。

#include <iostream>
#include<sys/time.h>
#include <Eigen/Dense>
 
using namespace std;
using namespace Eigen;
 
int main()

   Matrix3f A;
   Vector3f b;
   A << 1,2,3,  4,5,6,  7,8,10;
   b << 3, 3, 4;

    struct timeval start, end;
    gettimeofday(&start, NULL);

   cout << "Here is the matrix A:\\n" << A << endl;
   cout << "Here is the vector b:\\n" << b << endl;
   Vector3f x = A.householderQr().solve(b);
   cout << "The solution is:\\n" << x << endl;

   gettimeofday(&end, NULL);
   int timeuse = 1000000 * (end.tv_sec - start.tv_sec) + end.tv_usec - start.tv_usec;
   printf("线性方程组计算耗时:  %d  us\\n", timeuse);

输出：

Here is the matrix A:
 1  2  3
 4  5  6
 7  8 10
Here is the vector b:
3
3
4
The solution is:
      -2
0.999998
       1
线性方程组计算耗时:  480  us

2.2 ColPivHouseholderQR

ColPivHouseholderQR速度和精度位于3个分解方法中的中间状态，是一个很好的折中方法。

• 代码举例

#include <iostream>
#include<sys/time.h>
#include <Eigen/Dense>
 
using namespace std;
using namespace Eigen;
 
int main()

   Matrix3f A;
   Vector3f b;
   A << 1,2,3,  4,5,6,  7,8,10;
   b << 3, 3, 4;

    struct timeval start, end;
    gettimeofday(&start, NULL);

   cout << "Here is the matrix A:\\n" << A << endl;
   cout << "Here is the vector b:\\n" << b << endl;
   Vector3f x = A.colPivHouseholderQr().solve(b);
   cout << "The solution is:\\n" << x << endl;

   gettimeofday(&end, NULL);
   int timeuse = 1000000 * (end.tv_sec - start.tv_sec) + end.tv_usec - start.tv_usec;
   printf("线性方程组计算耗时:  %d  us\\n", timeuse);

输出：

Here is the matrix A:
 1  2  3
 4  5  6
 7  8 10
Here is the vector b:
3
3
4
The solution is:
      -2
0.999997
       1
线性方程组计算耗时:  482  us

2.3 FullPivHouseholderQR

FullPivHouseholderQR分解是3个QR分解中精度最高的一种，但是它的速度也是最慢的。

#include <iostream>
#include<sys/time.h>
#include <Eigen/Dense>
 
using namespace std;
using namespace Eigen;
 
int main()

   Matrix3f A;
   Vector3f b;
   A << 1,2,3,  4,5,6,  7,8,10;
   b << 3, 3, 4;

    struct timeval start, end;
    gettimeofday(&start, NULL);

   cout << "Here is the matrix A:\\n" << A << endl;
   cout << "Here is the vector b:\\n" << b << endl;
   Vector3f x = A.fullPivHouseholderQr().solve(b);
   cout << "The solution is:\\n" << x << endl;

   gettimeofday(&end, NULL);
   int timeuse = 1000000 * (end.tv_sec - start.tv_sec) + end.tv_usec - start.tv_usec;
   printf("线性方程组计算耗时:  %d  us\\n", timeuse);

以上是关于线性方程理论说明和Eigen解线性方程求解方法汇总的主要内容，如果未能解决你的问题，请参考以下文章
 
 使用 Eigen 求解线性方程组
 现代控制理论课程实验一：线性系统状态空间分析与运动解
 如何用matlab求解齐次线性方程组？举个例子说明
 怎么去解多元一次方程组快
 Eigen解线性方程组
 我可以使用 Eigen 求解形式为 xA=b 的线性方程组，其中 A 是稀疏的吗？