C++库Eigen学习2:普通线性问题和分解

前情提要

线性代数和分解

求解线性方程

根据系数矩阵的性质,对速度和精度的要求,可以选择不同的分解方式

下面这个例子总是可以解出来,并且在速度和精度上面有所权衡

#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;
}
分解方式函数对系数矩阵的要求速度(小矩阵)速度(大矩阵)精度
PartialPivLUpartialPivLu()Invertible+++++
FullPivLUfullPivLu()None– –+++
HouseholderQRhouseholderQr()None+++++
ColPivHouseholderQRcolPivHouseholderQr()None++++
FullPivHouseholderQRfullPivHouseholderQr()None– –+++
CompleteOrthogonalDecompositioncompleteOrthogonalDecomposition()None++++
LLTllt()Positive definite+++++++
LDLTldlt()Positive or negative
semidefinite
++++++
BDCSVDbdcSvd()None+++
JacobiSVDjacobiSvd()None– – –+++

不同分解方式的比较

更完成的比较参见:https://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html

各方法的求解速度比较:https://eigen.tuxfamily.org/dox/group__DenseDecompositionBenchmark.html

非对称满秩的矩阵可以用PartialPivLU,对称且正定的矩阵可以用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;
}

最小二乘解

解欠约束、过约束的线性方程就需要最小二乘解,常用的方法是奇异值分解(SVD)

Eigen自带两个奇异值分解的求解器,推荐使用的是BDCSVD。因为它对大矩阵的适应性更好,且对小矩阵来说等效于JacobiSVD

#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;
}

出现一个问题:节数超过对象文件格式限制: 请使用 /bigobj 进行编译

解决方法:右键项目-属性-C/C++-命令行-其他选项里面加上/bigobj

求解这类问题的另一个方法是CompleteOrthogonalDecomposition,通常更快且差不多精确

如果满秩,选HouseHolderQR。如果满秩且良置(低条件数),用Cholesky decomposition(LLT)更好

检查矩阵是否奇异

检查误差是否在容许范围内

#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;
}

计算特征值和特征向量

自伴随矩阵用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;
}

计算矩阵的逆和行列式

逆常常可以用solve()更好地替代,而行列式并不是检查矩阵是否可逆的好方法,但是对于小矩阵(4*4),则不一定

虽然PartialPivLU和FullPivLU提供了inverse()和determinant()函数,但是应该直接对矩阵用这两个函数,这样算小矩阵(4*4)会避免分解,算的更快

#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;
}

将计算与构造分离

上述类的构造和计算可以分开进行,这样就不用每次计算都构造一个变量了

#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;
}

预分配内存

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

反映矩阵的秩的分解

这些分解方法会计算矩阵的秩

一般来说这些方法是算不满秩(奇异)矩阵最好的方法

这些分解方法提供了rank(),isInvertible()等函数

#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
}

这些方法可以自定义计算的精度阈值,更改阈值之后不需要重新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;
}

求解最小二乘系统

SVD分解

#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;
}

QR分解

需要满秩

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;

一般方法

精度受条件数影响很大

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;

就地矩阵分解

直接对输入矩阵进行修改,适用于特大矩阵、内存小的系统

使用输入矩阵来存放中间值

#include <iostream>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

int main()
{
MatrixXd A(2, 2);
A << 2, -1, 1, 3;
cout << "Here is the input matrix A before decomposition:\n" << A << endl;

PartialPivLU<Ref<MatrixXd> > lu(A);
cout << "Here is the input matrix A after decomposition:\n" << A << endl;

cout << "Here is the matrix storing the L and U factors:\n" << lu.matrixLU() << endl;

MatrixXd A0(2, 2); A0 << 2, -1, 1, 3;
VectorXd b(2); b << 1, 2;
VectorXd x = lu.solve(b);
cout << "Residual: " << (A0 * x - b).norm() << endl;

A << 3, 4, -2, 1;
x = lu.solve(b);
cout << "Residual: " << (A0 * x - b).norm() << endl;

A0 = A; // save A
lu.compute(A);
x = lu.solve(b);
cout << "Residual: " << (A0 * x - b).norm() << endl;

MatrixXd A1(2, 2);
A1 << 5, -2, 3, 4;
lu.compute(A1);
cout << "Here is the input matrix A1 after decomposition:\n" << A1 << endl;

x = lu.solve(b);
cout << "Residual: " << (A1 * x - b).norm() << endl;
}

支持就地分解的方法如下

  • class LLT
  • class LDLT
  • class PartialPivLU
  • class FullPivLU
  • class HouseholderQR
  • class ColPivHouseholderQR
  • class FullPivHouseholderQR
  • class CompleteOrthogonalDecomposition

发表评论

您的电子邮箱地址不会被公开。

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据