# Eigen学习笔记(12)-线性代数与矩阵分解

The problem：对于一般形式如下的线性系统：
[ Ax : = : b ]

The solution：解决上述方程的方式一般是将矩阵A进行分解，你可以根据A的形式选择不同的分解方法。当然最基本的方法是高斯消元法。

``````#include <iostream>
#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;
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;
}
``````

``````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
1
``````

Eigen内置的解线性方程组的算法如下表所示：

Decomposition Method Requirements on the matrix Speed (small-to-medium) Speed (large) Accuracy
PartialPivLU partialPivLu() Invertible ++ ++ +
FullPivLU fullPivLu() None – – +++
HouseholderQR householderQr() None ++ ++ +
ColPivHouseholderQR colPivHouseholderQr() None + +++
FullPivHouseholderQR fullPivHouseholderQr() None – – +++
CompleteOrthogonalDecomposition completeOrthogonalDecomposition() None + +++
LLT llt() Positive definite +++ +++ +
LDLT ldlt() Positive or negative semidefinite +++ + ++
BDCSVD bdcSvd() None +++
JacobiSVD jacobiSvd() None – – – +++

``````#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
Matrix2f A, b;
A << 2, -1, -1, 3;
b << 1, 2, 3, 1;
cout << "Here is the matrix A:\\n" << A << endl;
cout << "Here is the right hand side b:\\n" << b << endl;
Matrix2f x = A.ldlt().solve(b);
cout << "The solution is:\\n" << x << endl;
}
``````

``````Here is the matrix A:
2 -1
-1  3
Here is the right hand side b:
1 2
3 1
The solution is:
1.2 1.4
1.4 0.8
``````

# 2. 确认解是否真的存在

``````#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
MatrixXd A = MatrixXd::Random(100,100);
MatrixXd b = MatrixXd::Random(100,50);
MatrixXd x = A.fullPivLu().solve(b);
double relative_error = (A*x - b).norm() / b.norm(); // norm() is L2 norm
cout << "The relative error is:\\n" << relative_error << endl;
}
``````

``````The relative error is:
2.31495e-14
``````

# 3. 计算特征值和特征向量

SelfAdjointEigenSolver 可以通过`EigenSolver`或者`ComplexEigenSolver`方法适用于一般矩阵。

``````#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
Matrix2f A;
A << 1, 2, 2, 3;
cout << "Here is the matrix A:\\n" << A << endl;
if (eigensolver.info() != Success) abort();
cout << "The eigenvalues of A are:\\n" << eigensolver.eigenvalues() << endl;
cout << "Here's a matrix whose columns are eigenvectors of A \\n"
<< "corresponding to these eigenvalues:\\n"
<< eigensolver.eigenvectors() << endl;
}
``````

``````Here is the matrix A:
1 2
2 3
The eigenvalues of A are:
-0.236
4.24
Here's a matrix whose columns are eigenvectors of A
corresponding to these eigenvalues:
-0.851 -0.526
0.526 -0.851
``````

# 4. 计算逆矩阵和矩阵行列式

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

`PartialPivLU``FullPivLU`提供了`inverse()``determinant()`方法，你也可以直接对矩阵使用`inverse()``determinant()`方法。如果矩阵是一个很小的固定尺寸矩阵（至少是4*4），允许Eigen阻止使用LU分解，在这种小矩阵上使用公式将会更有效。

``````#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
Matrix3f A;
A << 1, 2, 1,
2, 1, 0,
-1, 1, 2;
cout << "Here is the matrix A:\\n" << A << endl;
cout << "The determinant of A is " << A.determinant() << endl;
cout << "The inverse of A is:\\n" << A.inverse() << endl;
}
``````

``````
Here is the matrix A:
1  2  1
2  1  0
-1  1  2
The determinant of A is -3
The inverse of A is:
-0.667      1  0.333
1.33     -1 -0.667
-1      1      1
``````

# 5. 最小二乘法求解

Eigen也提供了解最小二乘问题的解法，并给出两种实现，分别是`BDCSVD``JacobiSVD``BDCSVD`可以很好地处理大问题，而对于较小的问题则自动返回到`JacobiSVD`类。因此推荐使用的一种是`BDCSVD`。这两个类的`solve()`方法都可以解决最小二乘问题。

``````#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
MatrixXf A = MatrixXf::Random(3, 2);
cout << "Here is the matrix A:\\n" << A << endl;
VectorXf b = VectorXf::Random(3);
cout << "Here is the right hand side b:\\n" << b << endl;
cout << "The least-squares solution is:\\n"
<< A.bdcSvd(ComputeThinU | ComputeThinV).solve(b) << endl;
}
``````

``````Here is the matrix A:
0.68  0.597
-0.211  0.823
0.566 -0.605
Here is the right hand side b:
-0.33
0.536
-0.444
The least-squares solution is:
-0.67
0.314
``````

# 6. 将计算与构造分离

• 所有分解都有一个默认构造函数，
• 所有分解都有一个`compute（matrix）`方法来执行计算，并且可以在已计算的分解上再次调用该方法，重新初始化它。
``````#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
Matrix2f A, b;
LLT<Matrix2f> llt;
A << 2, -1, -1, 3;
b << 1, 2, 3, 1;
cout << "Here is the matrix A:\\n" << A << endl;
cout << "Here is the right hand side b:\\n" << b << endl;
cout << "Computing LLT decomposition..." << endl;
llt.compute(A);
cout << "The solution is:\\n" << llt.solve(b) << endl;
A(1,1)++;
cout << "The matrix A is now:\\n" << A << endl;
cout << "Computing LLT decomposition..." << endl;
llt.compute(A);
cout << "The solution is now:\\n" << llt.solve(b) << endl;
}
``````

``````Here is the matrix A:
2 -1
-1  3
Here is the right hand side b:
1 2
3 1
Computing LLT decomposition...
The solution is:
1.2 1.4
1.4 0.8
The matrix A is now:
2 -1
-1  4
Computing LLT decomposition...
The solution is now:
1  1.29
1 0.571
``````

``````HouseholderQR<MatrixXf> qr(50,50);
MatrixXf A = MatrixXf::Random(50,50);
qr.compute(A); // no dynamic memory allocation
``````

# 7. 揭示秩的分解

** Rank-revealing decompositions** 提供了`rank()`方法，`isInvertible()`，有一些还提供了计算kernel（null-space)和image(column-space)的方法，示例如下：

``````#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
Matrix3f A;
A << 1, 2, 5,
2, 1, 4,
3, 0, 3;
cout << "Here is the matrix A:\\n" << A << endl;
FullPivLU<Matrix3f> lu_decomp(A);
cout << "The rank of A is " << lu_decomp.rank() << endl;
cout << "Here is a matrix whose columns form a basis of the null-space of A:\\n"
<< lu_decomp.kernel() << endl;
cout << "Here is a matrix whose columns form a basis of the column-space of A:\\n"
<< lu_decomp.image(A) << endl; // yes, have to pass the original A
}
``````

``````Here is the matrix A:
1 2 5
2 1 4
3 0 3
The rank of A is 2
Here is a matrix whose columns form a basis of the null-space of A:
0.5
1
-0.5
Here is a matrix whose columns form a basis of the column-space of A:
5 1
4 2
3 3
``````

``````#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
Matrix2d A;
A << 2, 1,
2, 0.9999999999;
FullPivLU<Matrix2d> lu(A);
cout << "By default, the rank of A is found to be " << lu.rank() << endl;
lu.setThreshold(1e-5);
cout << "With threshold 1e-5, the rank of A is found to be " << lu.rank() << endl;
}
``````

``````By default, the rank of A is found to be 2
With threshold 1e-5, the rank of A is found to be 1
``````

### 参考：

“Eigen学习之简单线性方程与矩阵分解”