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

前情提要

C++库Eigen学习1：普通矩阵和数组操作

线性代数和分解

求解线性方程

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

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

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

分解方式	函数	对系数矩阵的要求	速度（小矩阵）	速度（大矩阵）	精度
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	–	– – –	+++

不同分解方式的比较

更完成的比较参见：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