线性方程组/矩阵方程求解(方法汇总)

文章目录

前言
1 线性方程组（矩阵方程）
2 QR分解
- 2.1 ColPivHouseholderQR分解
- 2.2 HouseholderQR分解
- 2.3 FullPivHouseholderQR分解
3 LLT分解
3 LDLT分解
4 LU分解
- 4.1 LU分解
- 4.2 PartialPivLu分解
- 4.3 FullPivLu分解
5 SVD分解
- 5.1 BDCSVD分解
- 5.2 JacobiSVD分解
6 公式法（x=A−1bx=A^{-1}bx=A−1b）

前言

Eigen 是一个基于C++模板的线性代数库，提供了很多求解矩阵方程的方式，如比如LU分解、QR分解、SVD分解等。

本文以实用性为主，重点讨论几种求解矩阵方程的方法，以及不同方法的优缺点。对于每种方法的原理不做过多介绍。

下面给出不同分解方式的对比图，主要是计算速度和精度的对比，同时不同的分解方式对系数矩阵的形式也有不同的要求。

1 线性方程组（矩阵方程）

线性方程组是各个方程关于未知量均为一次的方程组，可以采用矩阵方程的形式表示。比如已知一个线性方程组
{a1x1+b1x2+c1x3=d1a2x1+b2x2+c2x3=d2a3x1+b3x2+c3x3=d3\begin{cases} 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{cases}⎩⎪⎨⎪⎧a1x1+b1x2+c1x3=d1a2x1+b2x2+c2x3=d2a3x1+b3x2+c3x3=d3
写成矩阵形式
Ax=bAx=bAx=b
其中，

A=[a1b1c1a2b2c2a3b3c3]A=\begin{bmatrix}a_1&b_1&c_1 \\a_2&b_2&c_2\\a_3&b_3&c_3\\\end{bmatrix}A=⎣⎡a1a2a3b1b2b3c1c2c3⎦⎤，x=[x1x2x3]x=\begin{bmatrix}x_1\\x_2\\x_3\\\end{bmatrix}x=⎣⎡x1x2x3⎦⎤，d=[d1d2d3]d=\begin{bmatrix}d_1\\d_2\\d_3\\\end{bmatrix}d=⎣⎡d1d2d3⎦⎤

以线性方程组
{x1+2x2+3x3=34x1+5x2+6x3=37x1+8x2+10x3=4\begin{cases} x_1+2x_2+3x_3=3\\ 4x_1+5x_2+6x_3=3\\ 7x_1+8x_2+10x_3=4\\ \end{cases}⎩⎪⎨⎪⎧x1+2x2+3x3=34x1+5x2+6x3=37x1+8x2+10x3=4
为例，进行求解。写成矩阵形式为
[1234567810][x1x2x3]=[334]\begin{bmatrix}1&2&3 \\4&5&6\\7&8&10\\\end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\\\end{bmatrix}=\begin{bmatrix}3\\3\\4\\\end{bmatrix}⎣⎡1472583610⎦⎤⎣⎡x1x2x3⎦⎤=⎣⎡334⎦⎤

可以选择不同的分解方式，这取决于矩阵A的形式，也取决于求解的速度和精度。

2 QR分解

QR分解类中的solution()方法也可以用来求解矩阵方程，有三个QR分解类：

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

2.1 ColPivHouseholderQR分解

ColPivHouseholderQR分解是一个很好的折衷方案，因为它适用于所有矩阵，同时速度非常快。

代码：

#include <iostream>
#include <Eigen/Dense>using namespace std;
using namespace Eigen;int main()
{Matrix3d A;    //系数矩阵AVector3d b;  //常数矩阵bVector3d x;  //未知数矩阵xA << 1,2,3,  4,5,6,  7,8,10;b << 3, 3, 4;cout << "系数矩阵A:\n" << A << endl;cout << "常数矩阵b:\n" << b << endl;x = A.colPivHouseholderQr().solve(b);cout << "未知数矩阵x:\n" << x << endl;return 0;
}

输出结果：

系数矩阵A:1  2  34  5  67  8 10
常数矩阵b:
3
3
4
未知数矩阵x:
-211

2.2 HouseholderQR分解

相比于ColPivHouseholderQR分解，HouseholderQR分解计算速度更快，但精度不如ColPivHouseholderQR分解。

代码

#include <iostream>
#include <Eigen/Dense>using namespace std;
using namespace Eigen;int main()
{Matrix3d A;    //系数矩阵AVector3d b;  //常数矩阵bVector3d x;  //未知数矩阵xA << 1,2,3,  4,5,6,  7,8,10;b << 3, 3, 4;cout << "系数矩阵A:\n" << A << endl;cout << "常数矩阵b:\n" << b << endl;x = A.householderQr().solve(b);cout << "未知数矩阵x:\n" << x << endl;return 0;
}

输出结果：

系数矩阵A:1  2  34  5  67  8 10
常数矩阵b:
3
3
4
未知数矩阵x:
-211

2.3 FullPivHouseholderQR分解

FullPivHouseholderQR分解相比于前两种QR分解，速度更慢，但精度更高。

代码：

#include <iostream>
#include <Eigen/Dense>using namespace std;
using namespace Eigen;int main()
{Matrix3d A;    //系数矩阵AVector3d b;  //常数矩阵bVector3d x;  //未知数矩阵xA << 1,2,3,  4,5,6,  7,8,10;b << 3, 3, 4;cout << "系数矩阵A:\n" << A << endl;cout << "常数矩阵b:\n" << b << endl;x = A.fullPivHouseholderQr().solve(b);cout << "未知数矩阵x:\n" << x << endl;return 0;
}

输出结果：

系数矩阵A:1  2  34  5  67  8 10
常数矩阵b:
3
3
4
未知数矩阵x:
-211

3 LLT分解

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

代码：

#include <iostream>
#include <Eigen/Dense>using namespace std;
using namespace Eigen;int main()
{Matrix3d A;    //系数矩阵AVector3d b;  //常数矩阵bVector3d x;  //未知数矩阵xA << 1,2,3,  4,5,6,  7,8,10;b << 3, 3, 4;cout << "系数矩阵A:\n" << A << endl;cout << "常数矩阵b:\n" << b << endl;x = A.llt().solve(b);cout << "未知数矩阵x:\n" << x << endl;return 0;
}

输出结果：

系数矩阵A:1  2  34  5  67  8 10
常数矩阵b:
3
3
4
未知数矩阵x:4.4556
-0.3184-0.026

3 LDLT分解

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

代码：

#include <iostream>
#include <Eigen/Dense>using namespace std;
using namespace Eigen;int main()
{Matrix3d A;    //系数矩阵AVector3d b;  //常数矩阵bVector3d x;  //未知数矩阵xA << 1,2,3, 4,5,6,  7,8,10;b << 3, 3, 4;cout << "系数矩阵A:\n" << A << endl;cout << "常数矩阵b:\n" << b << endl;x = A.ldlt().solve(b);cout << "未知数矩阵x:\n" << x << endl;return 0;
}

输出结果：

系数矩阵A:1  2  34  5  67  8 10
常数矩阵b:
3
3
4
未知数矩阵x:
-0.2068970.379310.241379

4 LU分解

4.1 LU分解

代码：

#include <iostream>
#include <Eigen/Dense>using namespace std;
using namespace Eigen;int main()
{Matrix3d A;    //系数矩阵AVector3d b;  //常数矩阵bVector3d x;  //未知数矩阵xA << 1,2,3,  4,5,6,  7,8,10;b << 3, 3, 4;cout << "系数矩阵A:\n" << A << endl;cout << "常数矩阵b:\n" << b << endl;x = A.lu().solve(b);cout << "未知数矩阵x:\n" << x << endl;return 0;
}

输出结果：

系数矩阵A:1  2  34  5  67  8 10
常数矩阵b:
3
3
4
未知数矩阵x:
-211

4.2 PartialPivLu分解

代码：

#include <iostream>
#include <Eigen/Dense>using namespace std;
using namespace Eigen;int main()
{Matrix3d A;    //系数矩阵AVector3d b;  //常数矩阵bVector3d x;  //未知数矩阵xA << 1,2,3,  4,5,6,  7,8,10;b << 3, 3, 4;cout << "系数矩阵A:\n" << A << endl;cout << "常数矩阵b:\n" << b << endl;x = A.partialPivLu().solve(b);cout << "未知数矩阵x:\n" << x << endl;return 0;
}

输出结果：

系数矩阵A:1  2  34  5  67  8 10
常数矩阵b:
3
3
4
未知数矩阵x:
-211

4.3 FullPivLu分解

代码：

#include <iostream>
#include <Eigen/Dense>using namespace std;
using namespace Eigen;int main()
{Matrix3d A;    //系数矩阵AVector3d b;  //常数矩阵bVector3d x;  //未知数矩阵xA << 1,2,3,  4,5,6,  7,8,10;b << 3, 3, 4;cout << "系数矩阵A:\n" << A << endl;cout << "常数矩阵b:\n" << b << endl;x = A.fullPivLu().solve(b);cout << "未知数矩阵x:\n" << x << endl;return 0;
}

生成失败：

error C1128: 节数超过对象文件格式限制: 请使用 /bigobj 进行编译

5 SVD分解

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

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

5.1 BDCSVD分解

代码：

#include <iostream>
#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;cout << "系数矩阵A:\n" << A << endl;cout << "常数矩阵b:\n" << b << endl;x = A.bdcSvd(ComputeThinU | ComputeThinV).solve(b);cout << "未知数矩阵x:\n" << x << endl;return 0;
}

输出结果：

系数矩阵A:
1 2
4 5
7 8
常数矩阵b:
3
3
4
未知数矩阵x:-2.5
2.66667

5.2 JacobiSVD分解

代码：

#include <iostream>
#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;cout << "系数矩阵A:\n" << A << endl;cout << "常数矩阵b:\n" << b << endl;x = A.jacobiSvd(ComputeThinU | ComputeThinV).solve(b);cout << "未知数矩阵x:\n" << x << endl;return 0;
}

输出结果：

系数矩阵A:
1 2
4 5
7 8
常数矩阵b:
3
3
4
未知数矩阵x:-2.5
2.66667

6 公式法（x=A−1bx=A^{-1}bx=A−1b）

虽然逆（inverse）和行列式（determinant）是基本的数学概念，但在线性代数中，它们不像在纯数学中那么流行。逆运算通常被 solve 运算所代替，行列式通常不是检查矩阵是否可逆的好方法。

但是，对于非常小的矩阵，逆矩阵和行列式是非常有用的。

虽然某些分解（如PartialPivLU和FullPivLU）提供了inverse()和determinate()方法，但也可以直接对矩阵调用inverse()和determinate()。如果矩阵是一个非常小的（最多4x4）静态矩阵，Eigen可使用公式 x=A−1bx=A^{-1}bx=A−1b，比矩阵分解的方式更有效。

代码：

#include <iostream>
#include <Eigen/Dense>using namespace std;
using namespace Eigen;int main()
{Matrix3d A;    //系数矩阵AVector3d b;  //常数矩阵bVector3d x;  //未知数矩阵xA << 1,2,3,  4,5,6,  7,8,10;b << 3, 3, 4;cout << "系数矩阵A:\n" << A << endl;cout << "常数矩阵b:\n" << b << endl;if (A.determinant() != 0){//矩阵行列式大于0，x = A.inverse()*b;cout << "未知数矩阵x:\n" << x << endl;}else{cout << "矩阵行列式小于0，矩阵方程无解！\a\n" << endl;}return 0;
}

输出结果：

系数矩阵A:1  2  34  5  67  8 10
常数矩阵b:
3
3
4
未知数矩阵x:
-211