本文详细介绍了Eigen库中解决线性方程组的不同方法，包括QR分解的HouseholderQR、ColPivHouseholderQR和FullPivHouseholderQR，LLT分解、LDLT分解、LU分解的partialPivLu和fullPivLu，以及SVD分解的BDCSVD和JacobiSVD。通过示例代码展示了各种方法的计算速度和精度，为选择合适的解算方法提供了参考。
文章目录
1. 线性方程组(矩阵方程)理论
线性方程组是各个方程关于未知量均为一次的方程组，方程表示如下：
{
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
\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.
⎩⎨⎧ a1x1+b1x2+c1x3=d1 a2x1+b2x2+c2x3=d2 a3x1+b3x2+c3x3=d3
把线程方程写成矩阵形式如下：
A
x
=
b
Ax = b
Ax=b
其中，
A
=
(
a
1
b
1
c
1
a
2
b
2
c
2
a
3
b
3
c
3
)
,
x
=
(
x
1
x
2
x
3
)
,
d
=
(
d
1
d
2
d
3
)
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}
A=⎝⎛a1a2a3b1b2b3c1c2c3⎠⎞,x=⎝⎛x1x2x3⎠⎞,d=⎝⎛d1d2d3⎠⎞
Eigen中提供了丰富的线性方程求解方法，包括LU分解法，QR分解法，SVD（奇异值分解）、特征值分解等根据A的矩阵类型、对结果的精度要求以及计算速度的要求，可以选择不同的计算方式，下面一一进行介绍。
下图是Eigen一些分解方法的简介。
如果对矩阵很了解，那么可以很方便的选择一种合理的求解方法，比如如果矩阵A是满秩且非对称矩阵，那么可以选用PartialPivLU求解方法，如果知道你的矩阵是对称正定的，那么选用LLT或者LDLT分解是一个很好的选择。下图是一些求解方法的速度。
2. QR分解
QR（正交三角）分解法是求一般矩阵全部特征值的最有效并广泛应用的方法，一般矩阵先经过正交相似变化成为Hessenberg矩阵，然后再应用QR方法求特征值和特征向量。它是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R，所以称为QR分解法，与此正规正交矩阵的通用符号Q有关。
A
=
Q
R
A = QR
A=QR
其中Q是正交矩阵（或酉矩阵），即
Q
Q
T
=
1
QQ^T=1
QQT=1，
R
R
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
=
L
L
T
=
R
T
R
A=LL^T=R^TR
A=LLT=RTR
其中的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为：
A
=
U
D
V
T
=
(
Σ
0
0
0
)
A = UDV^T = \begin{pmatrix} \Sigma & 0 \\ 0 & 0 \end{pmatrix}
A=UDVT=(Σ000)
其中，
U
U
U为
m
×
m
m\times m
m×m的矩阵，
Σ
\Sigma
Σ为
m
×
n
m\times n
m×n的矩阵，
V
V
V为
n
×
n
n\times n
n×n的矩阵。
Σ
=
d
i
a
g
(
σ
1
,
σ
2
,
.
.
.
,
σ
r
)
\Sigma=diag({\sigma}_1,{\sigma}_2,...,{\sigma}_r )
Σ=diag(σ1,σ2,...,σr)，
σ
{\sigma}
σ为矩阵A的全部非零奇异值。
U
U
U和
V
V
V满足，
U
T
U
=
I
,
V
T
V
=
I
U^TU=I,V^TV = I
UTU=I,VTV=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.youkuaiyun.com/qq_41839222/article/details/96274251
扫一扫
1358
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
