Eigen学习笔记20：线性代数及其分解_selfadjointeigensolver

基本的线性求解问题：

你有一个方程组，你已经写成一个矩阵方程Ax=b 其中A和b是矩阵（特殊情况下b可以是向量）。

您想找到一个解x。

解决方案：You can choose between various decompositions,  取决于你的矩阵是什么一个样子，取决于你选择速度或准确性。

#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
0.999997
       1

colPivHouseholderQr()方法返回类ColPivHouseholderQR的对象。

由于此处的矩阵类型为Matrix3f，因此该行可能已替换为：

ColPivHouseholderQR<Matrix3f> dec(A);
Vector3f x = dec.solve(b);

To get an overview of the true relative speed of the different decompositions, check this benchmark .

所有这些分解都提供了一个solve()方法，该方法与上述示例中的方法一样。

例如，如果您的矩阵是正定的，则​​上表说明LLT或LDLT分解是一个很好的选择。

#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

For a much more complete table comparing all decompositions supported by Eigen (notice that Eigen supports many other decompositions), see our special page on this topic.

检查解决方案是否确实存在

Only you know what error margin you want to allow for a solution to be considered valid.

因此，Eigen允许您自己进行此计算，如以下示例所示：

#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;
      SelfAdjointEigenSolver<Matrix2f> eigensolver(A);
      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.236068
  4.23607
Here's a matrix whose columns are eigenvectors of A 
corresponding to these eigenvalues:
-0.850651 -0.525731
 0.525731 -0.850651

计算特征值和特征向量

You need an eigendecomposition here, see available such decompositions on this page. 

确保检查矩阵是否是self-adjoint,，就像在这些问题中经常发生的那样。

这是一个使用SelfAdjointEigenSolver的示例，可以使用EigenSolver或ComplexEigenSolver轻松地将其应用于一般矩阵。

特征值和特征向量的计算不一定会收敛，但是这种收敛失败的情况非常罕见。调用info()就是为了检查这种可能性。

#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;
   SelfAdjointEigenSolver<Matrix2f> eigensolver(A);
   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.236068
  4.23607
Here's a matrix whose columns are eigenvectors of A 
corresponding to these eigenvalues:
-0.850651 -0.525731
 0.525731 -0.850651

计算逆和行列式

首先，请确保您确实想要这个。

尽管逆和行列式是基本的数学概念，但在数值线性代数中，它们不如在纯数学中流行。

Inverse运算通常可以用solve()操作代替，而行列式通常不是检查矩阵是否可逆的好方法。

但是，对于非常 小的矩阵，上述条件是不正确的，并且逆和行列式可能非常有用。

尽管某些分解（例如PartialPivLU和FullPivLU）提供了inverse()和determinant()方法，但您也可以直接在矩阵上调用inverse()和determinant()。

如果矩阵的固定大小很小（最多为4x4），这将使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.666667         1  0.333333
  1.33333        -1 -0.666667
       -1         1         1

最小二乘求解

最小二乘求解的最准确方法是SVD分解。

Eigen提供了两种实现。推荐的对象是BDCSVD类，它可以很好地解决较大的问题，并自动退回到JacobiSVD类以解决较小的问题。对于这两个类，它们的resolve()方法都在进行最小二乘求解。

这是一个例子：

#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:
 1  2  1
 2  1  0
-1  1  2
The determinant of A is -3
The inverse of A is:
-0.666667         1  0.333333
  1.33333        -1 -0.666667
       -1         1         1

可能更快但可靠性较低的另一种方法是使用法线矩阵的Cholesky分解或QR分解。decomposition. Our page on least squares solving has more details.

将计算与构造分开

在以上示例中，在构造分解对象的同时计算了分解。但是，在某些情况下，您可能希望将这两件事分开，例如，如果在构造时不知道要分解的矩阵，则可能会需要将它们分开。或者如果您想重用现有的分解对象。

使之成为可能的原因是：

• 所有分解都有默认的构造函数，
• 所有分解都具有执行计算的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.28571
       1 0.571429

最后，您可以告诉分解构造函数预先分配存储空间以分解给定大小的矩阵，以便在随后分解此类矩阵时，不执行动态内存分配（当然，如果您使用的是固定大小的矩阵，则不存在动态内存分配发生）。只需将大小传递给分解构造函数即可完成，如下例所示：

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

Rank-revealing 分解

某些分解是Rank-revealing ，即能够计算矩阵的等级。这些通常也是在非满秩矩阵（在方形情况下表示奇异矩阵）的情况下表现最佳的分解。On this table you can see for all our decompositions whether they are rank-revealing or not.

Rank-revealing decompositions offer at least a rank() method. 

它们还可以提供方便的方法，例如isInvertible()，并且还提供一些方法来计算矩阵的内核（空空间）和图像（列空间），就像FullPivLU那样：

#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

当然，任何秩计算都取决于对任意阈值的选择，因为实际上没有浮点矩阵恰好是秩不足的。Eigen选择一个明智的默认阈值，该阈值取决于分解，但通常是对角线大小乘以机器ε。虽然这是我们可以选择的最佳默认值，但只有您知道您的应用程序的正确阈值是多少。您可以通过在调用rank()或需要使用此阈值的任何其他方法之前，在分解对象上调用setThreshold()来进行设置。分解本身（即compute()方法）与阈值无关。更改阈值后，无需重新计算分解。

#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