【Eigen】Chapter2 线性方程组求解与矩阵分解 Dense Linear Problems and Decompositions

花舞の街道

于 2022-09-23 16:42:40 发布

阅读量639

点赞数 1

分类专栏： Eigen 文章标签：矩阵线性代数算法

本文链接：https://blog.csdn.net/ASUNAchan/article/details/127013851

版权

Eigen 专栏收录该内容

4 篇文章 0 订阅

订阅专栏

（1）线性代数和分解

介绍如何求解线性系统，计算各种分解，例如LU（三角分解），QR，SVD（奇异值分解），本征分解…

1）基本线性求解
$A x = b$
该解决方案：可将各种分解之间进行选择，取决于你的矩阵一个样子，取决于你是否赞成速度或准确性。但是，让我们从一个适用于所有情况的示例开始，这是一个很好的折衷方案：

#include <iostream>
#include <Eigen/Dense>
 
int main()
{
   Eigen::Matrix3f A;
   Eigen::Vector3f b;
   A << 1,2,3,  4,5,6,  7,8,10;
   b << 3, 3, 4;
   std::cout << "Here is the matrix A:\n" << A << std::endl;
   std::cout << "Here is the vector b:\n" << b << std::endl;
   Eigen::Vector3f x = A.colPivHouseholderQr().solve(b);
   std::cout << "The solution is:\n" << x << std::endl;
}

在此示例中，colPivHouseholderQr()方法返回类ColPivHouseholderQR的对象。由于此处的矩阵类型为Matrix3f，因此该行可能已替换为：

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

在这里，ColPivHouseholderQR是具有选择列主元功能的QR分解。这是本教程的一个不错的折衷方案，因为它适用于所有矩阵，而且速度非常快。

这是一些其他分解表，您可以根据矩阵和要进行的权衡选择：

分解	方法	矩阵要求	速度（低到中）	速度（高）	准确性
PartialPivLU	partialPivLu()	可逆	++	++	+
FullPivLU	fullPivLu()	无	-	–	+++
HouseholderQr	householderQr()	无	++	++	+
ColPivHouseholderQr	colPivHouseholderQr()	无	+	-	+++
FullPivHouseholderQr	fullPivHouseholderQr()	无	-	–	+++
CompleteOrthogonalDecomposition	completeOrthogonalDecomposition()	无	+	-	+++
LLT	llt()	正定	+++	+++	+
LDLT	ldlt()	半正定/半负定	+++	+	++
BDCSVD	bdcSvd()	无	-	-	+++
JacobiSvd	jacobiSvd()	无	-	—	+++

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

例如，如果矩阵是正定的，则上表说明LLT或LDLT分解是一个很好的选择。这是一个示例，也说明使用通用矩阵（而非矢量）作为右手边是可能的。

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

2）最小二乘问题

在最小二乘意义上解决欠定或超定线性系统的最通用和最准确的方法是 SVD 分解。 Eigen 提供了两种实现。推荐的一个是 BDCSVD 类，它可以很好地解决大问题，并自动回退到 JacobiSVD 类来解决小问题。对于这两个类，他们的 solve() 方法在最小二乘意义上求解了线性系统。

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

CompleteOrthogonalDecomposition 是 SVD 的替代方法，它通常更快且准确度差不多。同样，如果对问题有更多了解，上表包含可能更快的方法。如果矩阵是满秩的，HouseHolderQR 是首选方法。如果矩阵是满秩且条件良好的，则对正规方程的矩阵使用 Cholesky 分解 (LLT) 会更快。

3）检查矩阵是否奇异

#include <iostream>
#include <Eigen/Dense>
 
using Eigen::MatrixXd;
 
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
   std::cout << "The relative error is:\n" << relative_error << std::endl;
}

4）计算特征值与特征向量

以下示例为使用SelfAdjointEigenSolver的示例，可以使用EigenSolver或ComplexEigenSolver轻松地将其应用于一般矩阵。特征值和特征向量的计算不一定会收敛，但是这种收敛失败的情况很少。调用info（）就是为了检查这种可能性。

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

5）计算方阵的逆和行列式

首先，确保你真的想要这个。虽然逆和行列式是基本的数学概念，但在数值线性代数中它们不如在纯数学中有用。逆计算通常由solve() 操作代替，而行列式通常不是检查矩阵是否可逆的好方法。但是，对于非常小的矩阵，上述情况可能不成立，逆矩阵和行列式可能非常有用。

虽然某些分解（例如 PartialPivLU 和 FullPivLU）提供了 inverse() 和 determinant() 方法，但也可以直接在矩阵上调用 inverse() 和 determinant()。如果矩阵具有非常小的固定大小（最多 4x4），这允许 Eigen 避免执行 LU 分解，而是使用对此类小矩阵更有效的公式。

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

6）秩

某些分解是秩揭示的，即能够计算矩阵的秩。这些通常也是在面对非满秩矩阵（在方形情况下表示奇异矩阵）时表现最佳的分解。在上述表上，可以看到我们所有的分解，无论它们是否显示排名。秩揭示分解至少提供了一个 rank() 方法。它们还可以提供方便的方法，例如 isInvertible()，有些还提供计算矩阵的核（零空间）和图像（列空间）的方法，如 FullPivLU 的情况：

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

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

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

（2）求解线性最小二乘系统

本节介绍如何使用本征求解线性最小二乘系统。一个超定方程组，例如Ax = b，没有解。在这种情况下，在差异Ax - b尽可能小的意义上，搜索最接近解的向量x是有意义的。该x称为最小二乘解（如果使用欧几里得范数）。本节讨论的三种方法是SVD分解，QR分解和正态方程。其中，SVD分解通常最准确但最慢，正则方程(normal equations)最快但最不准确，QR分解介于两者之间。

1）使用SVD分解

BDCSVD类中的solve()中方法可以直接用来求解线性最小二乘系统。仅计算奇异值（此类的默认值）是不够的。还需要奇异矢量，但是稀疏SVD分解足以计算最小二乘解：

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

2）使用QR分解

QR分解类中的resolve()方法还计算最小二乘解。

有3种QR分解类：HouseholderQR（无pivoting，因此快速但不稳定），ColPivHouseholderQR（列枢轴，因此较慢但更准确）和FullPivHouseholderQR（全枢轴，因此最慢且最稳定）。

MatrixXf A = MatrixXf::Random(3, 2);
VectorXf b = VectorXf::Random(3);
cout << "The solution using the QR decomposition is:\n"
     << A.colPivHouseholderQr().solve(b) << endl;

3）使用正规方程

找到Ax = b的最小二乘解等效于求解法线方程 $A^TAx = A^Tb$ 。如果矩阵A是病态的，那么这不是一个好方法，因为 $A^TA$ 的条件数是A的条件数的平方。这意味着与使用其他方法相比，使用正则方程式丢失的数字要多两倍。

病态矩阵。。。没学过

MatrixXf A = MatrixXf::Random(3, 2);
VectorXf b = VectorXf::Random(3);
cout << "The solution using normal equations is:\n"
     << (A.transpose() * A).ldlt().solve(A.transpose() * b) << endl;