关于Ax=b形式的矩阵分解与求解问题

对于一个线性方程\(Ax=b\)，其中A的维度为\(m\times n， m \geq n\)，\(x\)的维度为\(n\times 1\)，b的维度为\(m \times 1\)，则在A满秩的情况下，\(x\)有解且唯一。

对于高斯牛顿法，我们一般求解的是\(J^T J x = J^T b\)，其中\(J\)为残差\(b\)对状态\(x\)的雅克比矩阵，维度为\(m*n\)，这样的话其实更加简单了，因为线性方程的系数矩阵\(A\)变成了一个方阵，并且这个方阵是实对称且半正定的。

高斯消元法

高斯消元法的原理是，将矩阵A通过初等变换变成一个阶梯型的矩阵（变成一个\(n\times n\)的上三角矩阵），然后逐行求解；

QR分解

将矩阵A分解成一个正交矩阵和一个上三角矩阵的乘积\(A=QR\)，所以只需要求解\(Rx=Q^Tb\)即可，其中\(R\)为一个上三角矩阵；进行qr分解主要有三种方法：

• Househoder
  • 时间复杂度：\(\frac{2}{3}n^3\)
• Gram-Schmidt
  • 时间复杂度：\(\frac{3}{3}n^3\)
• Givens rotation
  • 时间复杂度：\(\frac{4}{3}n^3\)

Givens Rotation求解思路

另外一种消元思路，也是将A变换为一个上三角矩阵。这种思路是不断的对A进行二维的旋转，每次旋转会将一个位置的元素置为0，经过若干次之后将下三角部分所有的元素都变成了0，givens rotation也就达到目标了。

因为旋转矩阵是一种正交矩阵，因此givens rotation就是对原本的A不断的乘以旋转矩阵，最终将A转化成了一个上三角矩阵R，那么之前乘以的所有旋转矩阵累乘起来就是上面qr分解中的\(Q^T\)。

例子：

首先，对于一个二维矩阵块 \[ {\begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix}} \]

我们使用二维旋转矩阵 \[ {\begin{bmatrix} cos{\theta} & -sin{\theta} \\ sin{\theta} & cos{\theta} \end{bmatrix}} \]

左乘它，在满足\(\theta = atan(-\frac{a_{21}}{a_{11}})\)时，可以得到结果 \[ {\begin{bmatrix} b_{11} & b_{12} \\ 0 & b_{22} \end{bmatrix}} \] （满足其他角度条件可以将其他位置的元素变为0）。

在使用二维旋转矩阵不断的将下三角的元素变成0的过程中，需要遵循：

• 从左往右逐列进行；
• 每一列，从下往上执行；

对于一个维度\((m,n)\)的矩阵A，我们对其进行左乘一个正交矩阵Q应该是维度为\((m,m)\)的。因此假设我们要对A(m1,n1,2,2)这个矩阵块中的一个元素变为0，则Q矩阵应该首先被赋值成一个单位矩阵，然后将Q(m1,n1,2,2)矩阵块赋值为二维旋转矩阵（旋转矩阵由A(m1,n1,2,2)中要消去的变量确定）。按照上面的顺序操作可以使得一个被变为0的元素不会再变为非0值。

出现：1. openvins代码中MSCKF更新模块将测量误差方程投影到H_f的零空间上，消除路标的影响；2. openvins代码中MSCKF更新模块零空间投影完得到系统总的重投影测量误差方程后，通过givens rotation将H_x变为一个上三角方阵，res变为与其同维度向量，消减计算量。

Cholesky分解

将A分解为\(A = L L^T\)，其中L为下三角矩阵。然后，首先求解\(Lc = b\)，c为中间结果;再求解\(L^T x = c\)，得到最终结果。此过程要求A实对称且正定。

由于Cholesky分解对A的要求，所以一般情况下，对于方程\(Ax=b\)，会首先变成\((A^T A + \lambda I) x = A^T b\)再使用（刚好是LM优化的形式，因此LM优化时总会使用llt分解求解）。

LU分解

LU分解将A分解为一个下三角矩阵与一个上三角矩阵的乘积，如A= LU。然后可以通过两次高斯消元法进行求解。

事实上，并不是每个矩阵都有LU分解。而这种情况可以通过对A各行顺序重排进行解决，此时的分解叫做PLU分解，分解形式为A=PLU，其中P为置换矩阵（主要用来对行序进行交换的），L为下三角矩阵，U为上三角矩阵，此时其实可以这样理解，\(P^{-1} A = LU\)，即对A交换行序后再进行LU分解。

如果不止交换行序，还交换列序，则是A = PLUQ，Q是一个排列矩阵。