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

最新推荐文章于 2024-01-30 21:24:07 发布

非晚非晚

最新推荐文章于 2024-01-30 21:24:07 发布

阅读量3.3k

点赞数 5

分类专栏：数学理论与代码实现数学理论与工具文章标签：线性代数 eigen LU分解法 QR分解法 SVD分解法

本文链接：https://blog.csdn.net/QLeelq/article/details/122577589

版权

数学理论与工具同时被 2 个专栏收录

10 篇文章 27 订阅

订阅专栏

数学理论与代码实现

3 篇文章 0 订阅

订阅专栏

文章目录

1. 线性方程组(矩阵方程)理论
2. QR分解
3. LLT分解
4. LDLT分解
5. LU分解
6. SVD分解
- 6.1 BDCSVD分解
- 6.2 JacobiSVD分解
7. 逆矩阵法

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

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

$\left\{\begin{matrix} \ a_1x_1+b_1x_2+c_1x_3=d_1\\ \ a_2x_1+b_2x_2+c_2x_3=d_2\\ \ a_3x_1+b_3x_2+c_3x_3=d_3 \end{matrix}\right.$
把线程方程写成矩阵形式如下：
$A x = b$
其中，
$A=\begin{pmatrix} a_{1}& b_{1}& c_{1}\\ a_{2}& b_{2}& c_{2}\\ a_{3}& b_{3}& c_{3} \end{pmatrix},x =\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix},d=\begin{pmatrix}d_1\\d_2\\d_3\end{pmatrix}$

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

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

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

2. QR分解

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

$A = Q R$
其中Q是正交矩阵（或酉矩阵），即 $QQ^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);
}

输出：

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
       1
0.999999
线性方程组计算耗时:  131  us

3. LLT分解

LLT分解要求系数矩阵A是正定矩阵(Positive definite matrix)，是所有分解方式中速度最快的一种。用于非正定矩阵的分解时，难以得到正确的结果。

LLT分解即矩阵的Cholesky分解，又被称为平方根分解，是LDLT分解的一种特殊形式，即其中的D为单位矩阵。对称正定矩阵A可以分解成一个下三角矩阵L和L的转置LT相乘的形式：
$A=LL^T=R^TR$
其中的L是下三角矩阵，R是上三角矩阵。

#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.llt().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:
 4.4556
-0.3184
 -0.026
线性方程组计算耗时:  97  us

4. LDLT分解

LDLT分解要求系数矩阵A是正定矩阵(positive definite matrix)或者半负定矩阵(negative semi-definite matrix)。用于其他矩阵的分解时，难以得到正确的结果。

A为对称矩阵，且任意一K阶主子阵均不为0时，A有如下唯一的分解形式：

即L为下三角单位矩阵，D为对角矩阵。LDLT方法实际上是Cholesky分解法的改进（LLT分解需要开平方），用于求解线性方程组。

#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.ldlt().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:
 4.4556
-0.3184
 -0.026
线性方程组计算耗时:  100  us

5. LU分解

LU分解(LU Decomposition)是矩阵分解中最普通的一种，也是最经典的一种，它可以将一个矩阵分解为一个单位下三角矩阵和一个上三角矩阵的乘积。将所给的系数矩阵A转变成等价两个矩阵L和U的乘积，其中L和U分别是单位下三角矩阵和上三角矩阵。

5.1 LU分解

#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.lu().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
       1
0.999999
线性方程组计算耗时:  200  us

5.2 partialPivLu

#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.partialPivLu().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
       1
0.999999
线性方程组计算耗时:  390  us

5.3 fullPivLu

#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.fullPivLu().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
线性方程组计算耗时:  397  us

6. SVD分解

奇异值分解(Singular Value Decomposition，以下简称SVD)是在机器学习领域广泛应用的算法，它不光可以用于降维算法中的特征分解，还可以用于推荐系统，以及自然语言处理等领域。是很多机器学习算法的基石。SVD也是对矩阵进行分解，但是和特征分解不同，SVD并不要求要分解的矩阵为方阵。假设我们的矩阵A是一个m×n的矩阵，那么我们定义矩阵A的SVD为：
$UDV^T = \begin{pmatrix} \Sigma & 0 \\ 0 & 0 \end{pmatrix}$
其中， $U$ 为 $m\times m$ 的矩阵， $\Sigma$ 为 $m\times n$ 的矩阵， $V$ 为 $n\times n$ 的矩阵。 $\Sigma=diag({\sigma}_1,{\sigma}_2,...,{\sigma}_r )$ ， ${\sigma}$ 为矩阵A的全部非零奇异值。 $U$ 和 $V$ 满足， $U^TU=I,V^TV = I$

当方程个数大于未知数个数时，需要进行最小二乘求解，最有效的办法就是SVD分解，Eigen提供了两种分解方法，分别是BDCSVD和JacobiSVD。一般情况下会推荐使用BDCSVD，它可以很好地扩展大型矩阵，而对于较小的矩阵，它会自动退回到JacobiSVD类。

注意：采用SVD分解时，要用动态矩阵定义系数矩阵A、常数矩阵b、未知数矩阵x。

6.1 BDCSVD分解

#include <iostream>
#include<sys/time.h>
#include <Eigen/Dense>
 
using namespace std;
using namespace Eigen;
 
int main()
{
	MatrixXd A(3,2);	//系数矩阵A，3×2阶
	VectorXd b(3,1);	//常数矩阵b，3×1阶
	VectorXd x(3,1);	//未知数矩阵x，3×1阶

   A << 1, 2, 4, 5, 7, 8;
   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;
   x = A.bdcSvd(ComputeThinU | ComputeThinV).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
4 5
7 8
Here is the vector b:
3
3
4
The solution is:
   -2.5
2.66667
线性方程组计算耗时:  675  us

6.2 JacobiSVD分解

#include <iostream>
#include<sys/time.h>
#include <Eigen/Dense>
 
using namespace std;
using namespace Eigen;
 
int main()
{
	MatrixXd A(3,2);	//系数矩阵A，3×2阶
	VectorXd b(3,1);	//常数矩阵b，3×1阶
	VectorXd x(3,1);	//未知数矩阵x，3×1阶
	A << 1, 2,  4, 5,  7, 8;
	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;
   x = A.jacobiSvd(ComputeThinU | ComputeThinV).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
4 5
7 8
Here is the vector b:
3
3
4
The solution is:
   -2.5
2.66667
线性方程组计算耗时:  578  us

7. 逆矩阵法

Eigen 也提供了求逆矩阵和求矩阵行列式的算法，但是这两种算法对于大型矩阵来说都是非常不经济的算法，当需要对大型矩阵做这种的操作时，需要自己判断到底需不需这样做。但是对于小型矩阵则可以没有顾虑地使用。

#include <iostream>
#include<sys/time.h>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

int main()
{
    Matrix3d A; //系数矩阵A
    Vector3d b; //常数矩阵b
    Vector3d x; //未知数矩阵x
    A << 1, 2, 3, 4, 5, 6, 7, 8, 10;
    b << 3, 3, 4;
    cout << "系数矩阵A:\n"
         << A << endl;
    cout << "常数矩阵b:\n"
         << b << endl;

    struct timeval start, end;
    gettimeofday(&start, NULL);
    if (A.determinant() != 0)
    {
        //矩阵行列式大于0，
        x = A.inverse() * b;
        cout << "未知数矩阵x:\n"
             << x << endl;
    }
    else
    {
        cout << "矩阵行列式小于0，矩阵方程无解！\a\n"
             << endl;
    }

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

输出：

系数矩阵A:
 1  2  3
 4  5  6
 7  8 10
常数矩阵b:
3
3
4
未知数矩阵x:
-2
 1
 1
线性方程组计算耗时:  63  us

参考文档：
https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html
https://blog.csdn.net/qq_41839222/article/details/96274251

非晚非晚

关注

5
点赞
踩
40

收藏

觉得还不错? 一键收藏
打赏
0
评论
线性方程理论说明和Eigen解线性方程求解方法汇总

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

扫一扫