Eigen 矩阵的LDLT分解求解线性方程组

1.LDLT分解原理

  利用矩阵$A$的$LD{L}^T$分解来求解方程组$Ax=b$的方法称为$LD{L}^T$分解法。若对称矩阵$A$的各阶顺序主子式不为零则$A$可以唯一分解为$A=LD{L}^T$。其中$L$和$D$的形式如下
$$L=\left|\begin{array}{ccccc}1& & & & \\ l_{21}& 1& & & \\ l_{31}& l_{32}& 1& & \\ \cdots & \cdots & \cdots & \ddots & \\ l_{n1}& l_{n2}& \cdots & l_{nn-1}& 1\end{array}\right|$$

$$D=\left|\begin{array}{cccc}{d}_1& & & \\ & d_2& & \\ & & \ddots & \\ & & & d_n\end{array}\right|$$
  当进行对称线性方程组$Ax=b$的求解时由$A=LD{L}^T$可得$(LD{L}^T)x=b$即$LD{L}^{T}x=b$。我们分别令
$$\begin{gathered} z=L^{T}x \\ y=Dz \end{gathered}$$
则有
$$Ly=b$$
因此可分步解得$x_i$。
$$\begin{gathered} y_i=b_i-\sum _{k=1}^{i-1}{l}_{ik}{y}_k,i=1,2,\cdots ,n \\ {z}_i=y_i/d_i,i=n,n-1,\cdots ,1 \\ {x}_i=z_i-\sum _{k=i+1}^n{l}_{ki}{x}_k,i=n,n-1,\cdots ,1 \end{gathered}$$

2.Eigen库实现

2.1问题定义

  用$LD{L}^T$法求解线性方程组$Ax=b$其中$A$和$b$如下
$$A=\left|\begin{array}{ccc}1& 1/2& 1/2\\ 1/2& 1& 1/2\\ 1/2& 1/2& 1\end{array}\right|$$

$$b=\left|\begin{array}{c}1\\ -2\\ 3\end{array}\right|$$

2.2代码实现

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
	Matrix3d A;//3x3的矩阵
	Vector3d b;//3x1的列向量
	A << 1, 0.5, 0.5, 0.5, 1, 0.5, 0.5, 0.5, 1;//给A赋值
	b << 1, -2, 3;//给b赋值
	Vector3d x = A.ldlt().solve(b);//求解Ax=b
	cout <<"x^T=[" <<x.transpose()<<"]" << endl;
	return 0;
}

2.3输出结果

x^T=[ 1 -5  5]