Eigen3学习笔记——密集线性问题和分解

密集线性问题和分解

线性代数及其分解

基本线性方程求解

问题的数学描述，$A,b$为矩阵（$b$当然可以是向量），$x$是解：
$$Ax=b$$
EIgen 解决着这种线性方程后很多方法，但调用的格式是一样的，以其中一个为例：

Matrix3f A;
Vector3f b;
A << 1,2,3,  4,5,6,  7,8,10;
b << 3, 3, 4;
//注意调用格式
Vector3f x = A.colPivHouseholderQr().solve(b);
//A.colPivHouseholderQr()返回的是一个ColPivHouseholderQR类对象
//因此，也可以这么搞
ColPivHouseholderQR<Matrix3f> dec(A);
Vector3f x = dec.solve(b);

以上是统一的调用格式，具体有哪些方法可以进行线性方程组求解，见下表：

名称（不一定对）	分解的方法类（Eigen中真实存在的类）	函数	对于矩阵的要求
部分主元LU分解	PartialPivLU	partialPivLu()	可逆
全主元LU分解	FullPivLU	fullPivLu()	无
Householder QR分解	HouseholderQR	householderQr()	无
列主元Householder QR分解	ColPivHouseholderQR	colPivHouseholderQr()	无
全主元Householder QR分解	FullPivHouseholderQR	fullPivHouseholderQr()	无
完全正交分解COD	CompleteOrthogonalDecomposition	completeOrthogonalDecomposition()	无
标准Cholesky分解	LLT	llt()	正定
鲁棒Cholesky分解	LDLT	ldlt()	正半定或负半定
双对角分治SVD分解	BDCSVD	bdcSvd()	无
雅各比svd分解	JacobiSVD	jacobiSvd()	无

最小二乘求解

求解最小二乘最通用和最准确的方法是 SVD 分解，

1. 推荐使用的是BDCSVD，
2. 另一种替代的方法是CompleteOrthogonalDecomposition
3. 也可以使用QR分解的HouseholderQR，ColPivHouseholderQR，FullPivHouseholderQR
4. 还可以使用正规方程$A^{T}Ax=A^{T}b$：(A.transpose()*A).ldlt().solve(A.transpose()*b)

计算特征值和特征向量

1. SelfAdjointEigenSolver
2. ComplexEigenSolver
3. EigenSolver
4. GeneralizedSelfAdjointEigenSolver

计算逆和行列式

虽然取逆和行列式是数学上的基本概念，但是对于大矩阵在实现上不不可取，当然Eigen 依然提供这种方法，针对的是小矩阵。

1. A.determinant()取行列式
2. A.inverse()取逆

将分解与构造对象分离

很多情况下，并不是像上面那样，直接对矩阵进行构造然后求解，而是分开进行，比如：你在分解A的时候不知道b是谁。为此，EIgen 是很灵活的，每一个分解器都有默认的构造函数，所以我们可以先创建，在分别提供所要分解的矩阵A和另一个矩阵b。

Matrix2f A, b;
LLT<Matrix2f> llt;		//注意此处求解器构造的类型
A << 2, -1, -1, 3;
b << 1, 2, 3, 1;
//先A后b
llt.compute(A);
llt.solve(b);
//并且求解器可以通用
A(1,1)++;
llt.compute(A);
llt.solve(b);

//对于动态大小的求解器也可以预先分配内存
HouseholderQR<MatrixXf> qr(50,50);
MatrixXf A = MatrixXf::Random(50,50);
qr.compute(A); // no dynamic memory allocation

计算矩阵的秩

有许多分解本身即使通过矩阵的秩进行的，因此可以在那些分解的基础上使用.rank()求得矩阵的秩。

Matrix3f A;
A << 1, 2, 5,
	2, 1, 4,
	3, 0, 3;
	
FullPivLU<Matrix3f> lu_decomp(A);
//利用求解器计算矩阵的秩
lu_decomp.rank();
//利用求解器计算矩阵的零空间
lu_decomp.kernel();
//利用求解器计算矩阵的列空间
lu_decomp.image(A);
//因为计算是有精度的，求解器的阈值大小将影响结果
lu.setThreshold(1e-5);

注意：Eigen 中的秩的求解与一般的数学思路是不一样，数学上可能直接就用矩阵进行计算从A直接得到秩，然而Eigen 中是先分解后给出秩。我觉得很奇特。

矩阵分解的内置

Eiegn 3.3以后对于一些分解提供了矩阵分解内置的特性，上面所用到的求解器在得到A后，A依然存在，而求解器记录了分解的矩阵，相当于开辟了新空间。但是当采用矩阵分解内置后，求解器将分解的矩阵内置到A中，A被破坏，这样节省空间，对于嵌入式系统很有用。使用这种方法需要在矩阵类型实例化的时候使用Ref<>。

MatrixXd A(2,2);
A << 2, -1, 1, 3;
PartialPivLU<Ref<MatrixXd> > lu(A);			
//矩阵分解内置，此后A被破坏，lu中记录了分解结果
//以后lu求解器的分解结果都会放在这个内存位置

MatrixXd A1(2,2);
A1 << 5,-2,3,4;
lu.compute(A1);		//利用lu取分解A1，A1不会被破坏，分解的值依然在之前A的内存位置。

支持矩阵分解内置的类

LLT，LDLT，PartialPivLU，FullPivLU，HouseholderQR，ColPivHouseholderQR，FullPivHouseholderQR，CompleteOrthogonalDecomposition.

官方给出的具体表格

1. 矩阵分解表（Catalogue of dense decompositions）
2. 矩阵分解效率基准表（Benchmark of dense decompositions）