加速LU分解和非奇异矩阵求逆

1. 背景介绍

上(下)三角矩阵有许多性质使得计算机可以方便地对它们进行各种线性代数基础运算。例如对三角矩阵的求逆非常简单易理解。因此在计算非奇异矩阵的逆时常将原矩阵$A$分解为一个下三角矩阵$L$（$L$矩阵的对角线元素为$1$，对角线以上部分元素为$0$）和一个上三角矩阵$U$（$U$矩阵对角线以下部分元素为$0$）相乘的形式，再分别求下三角矩阵的逆$L^{-1}$和上三角矩阵的逆$U^{-1}$，再将它们相乘获得最终的结果。

理论上利用LU分解求矩阵逆分为三步:

• 对矩阵进行LU分解：$A=LU$；
• 求解线性方程组：$LX=I$，$X=L^{-1}$；
• 求解线性方程组：$UB=L^{-1}$, $B=U^{-1}{L}^{-1}=A^{-1}$；

LU分解可以用于任意宽高的矩阵，本文只讨论方阵。

2. 理论基础

在进行$LU$分解时，需要构造一个下三角矩阵和一个上三角矩阵，这时需要满秩变换来逐向量地将原矩阵变换为上三角矩阵$U$，并记录变换矩阵$L$。

2.1 高斯消元

2.1节介绍$LU$分解的关键步骤-高斯消元。假设此时需要将$n$维矩阵$A$的第$k$($1\le k\le n$)列向量$a_{\cdot ,k}={\left[a_{1,k},a_{2,k},\cdots ,a_{n,k}\right]}^T$中的第$k+1$个元素到第$n$个元素变换为$0$，需要构造变换矩阵$M_k$。首先需要构造高斯消元向量${\tau }_k$：
$${\tau }_k^T=\left[0,\cdots ,0,{\tau }_{k+1},\cdots ,{\tau }_n\right],其中{\tau }_i=\frac{a_{i,k}}{a_{k,k}},i=k+1,\cdots ,n$$
再定义辅助矩阵$E_k$：
$$E_k=\left[\begin{array}{ccccc}0& \cdots & 1& \cdots & 0\\ ⋮& & ⋮& & ⋮\\ 0& \cdots & 1& \cdots & 0\\ ⋮& & ⋮& & ⋮\\ 0& \cdots & 1& \cdots & 0\end{array}\right],E_k中第k列全为1,其他元素为0$$
于是可以得到高斯消元矩阵$M_k$：
$$M_k=I-{\tau }_k\odot E_k,(\odot 表示点乘)$$
最后对$a_{\cdot ,k}$进行高斯消元：
$$M_k\cdot a_{\cdot ,k}=\left[\begin{array}{cccccc}1& \cdots & 0& 0& \cdots & 0\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ 0& \cdots & 1& 0& \cdots & 0\\ 0& \cdots & -{\tau }_{k+1}& 1& \cdots & 0\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ 0& \cdots & -{\tau }_n& 0& \cdots & 1\end{array}\right]\left[\begin{array}{c}{a}_{1,k}\\ ⋮\\ a_{k,k}\\ a_{k+1,k}\\ ⋮\\ a_{n,k}\end{array}\right]=\left[\begin{array}{c}{a}_{1,k}\\ ⋮\\ a_{k,k}\\ 0\\ ⋮\\ 0\end{array}\right]\text{\hspace{1em}(2.1.1)}$$

2.2 矩阵三角化

对$n$维方阵$A$的每一列从左到右重复$n$次高斯消元便可以得到上三角矩阵$U$:
$$M_n(M_{n-1}(\cdots (M_{1}A))=U\text{\hspace{1em}(2.1.1)}$$
$L$矩阵即是$n$个高斯消元矩阵的乘积:
$$A=(M_n{M}_{n-1}\cdots M_1{)}^{-1}U=LU\text{\hspace{1em}(2.1.2)}$$

$$L=M_1^{-1}{M}_2^{-1}\cdots M_n^{-1}\text{\hspace{1em}(2.1.3)}$$

2.3 LU分解的一般形式

设非奇异矩阵$A\in {\mathbb{R}}^{n\times n}$
$$A=\left[\begin{array}{cc}\alpha & {\omega }^T\\ v& B\end{array}\right],其中\alpha \in \mathbb{R},{\omega }^T\in {\mathbb{R}}^{1\times (n-1)},v\in {\mathbb{R}}^{(n-1)\times 1},B\in {\mathbb{R}}^{(n-1)\times (n-1)}$$
那么高斯消元的第$1$步分解结果是：
$$\left[\begin{array}{cc}{\alpha }_1& {\omega }_1^T\\ v_1& B_1\end{array}\right]=\left[\begin{array}{cc}1& 0\\ z_1/{\alpha }_1& I_1\end{array}\right]\left[\begin{array}{cc}1& 0\\ 0& B_1-z_1{w}_1^T/{\alpha }_1\end{array}\right]\left[\begin{array}{cc}{\alpha }_1& {\omega }_1^T\\ 0& I_1\end{array}\right]\text{\hspace{1em}(2.3.1)}$$
第$k,k=2,3,\cdots ,n$步可以重复以上步骤分解$B_{k-1}-z_{k-1}{w}_{k-1}^T/{\alpha }_{k-1}=L_{k-1}{U}_{k-1}$。那么：
$$L_{k-1}{U}_{k-1}=\left[\begin{array}{cc}1& 0\\ z_k/{\alpha }_k& I_k\end{array}\right]\left[\begin{array}{cc}1& 0\\ 0& B_k-z_k{w}_k^T/{\alpha }_k\end{array}\right]\left[\begin{array}{cc}{\alpha }_k& {\omega }_k^T\\ 0& I_k\end{array}\right]=\left[\begin{array}{cc}1& 0\\ z_k/{\alpha }_k& L_k\end{array}\right]\left[\begin{array}{cc}{\alpha }_k& {\omega }_k^T\\ 0& U_k\end{array}\right]\text{\hspace{1em}(2.3.2)}$$

3. 算法实现

对于方阵$A$，假设在已经过$k-1$步高斯消元，那么下一步是对$A(:,k)$进行高斯消元。再解下来便是利用高斯消元继续更新$A(k+1:n,k+1:n)$。由此可以推导出3.1节中的算法。

3.1 朴素的LU分解算法

假设$A\in {\mathbb{R}}^{n\times n}$非奇异，$L$矩阵是下三角矩阵（对角线元素为$1$），$U$是上三角矩阵，在进行$LU$分解之后，$A$的下三角部分被$L$覆盖（不包括对角线），上三角部分被$U$覆盖（包括对角线）。以下是算法伪代码：

  for $k=1:n-1$

​    $\rho =k+1:n$

​    $A(\rho ,k)=A(\rho ,k)/A(k,k)$

​    $A(\rho ,\rho )=A(\rho ,\rho )-A(\rho ,k)\cdot A(k,\rho )$

​  end

朴素的$LU$分解算法在C语言实现时有三层循环，那么排列组合$A_3^3$之后一共有$6$个版本。假设矩阵在内存中以数组的形式按行顺序存储，那么一种版本的C语言代码为：

void LU_decomposition(double* A, const size_t n)
{
    size_t i, j, k;
    double r;
    for (k=0; k<n-1; k++)
    {
        for (i=k+1; i<n; i++)
        {
            A[i * n + k] =  A[i * n + k] / A[k * n + k];
            for (j=k+1; j<n; j++)
                A[i * n + j] = A[i * n + j] - A[i * n + k] * A[k * n + j];
        }
    }
}

注意此处编译器优化：
例如在第9行的C语言代码中，会取出A元素的第i * n + k个元素，编译器判断每一次取数的地址是等间隔递增的，因此在程序编译完成后的汇编代码会使用简单的加法替代乘法。在函数上下文中，参数A和n的值在寄存器r0和r1中，假设A[i * n + k]元素的地址已在r6寄存器中：

MUL r5, r1, #8（行间隔字节量，双精度浮点数占8字节）

…

ADD r6, r6, r5

FLDD d0, [r6]

…

3.2 Gaxpy vs 外积

在工程应用中，两种计算量相同的算法常常会有不同的数据迁移量。考虑一下$n\times n$的Gaxpy计算和外积计算：
$$y=y+Ax,$$

$$A=A+y{x}^T,A\in {\mathbb{R}}^{n\times n},x,y\in {\mathbb{R}}^{n\times 1}$$

比较两种计算量相同算法的计算形式，用$\leftarrow$表示内存和寄存器之间的读写，对于Gaxpy：

​  $r_x\leftarrow x,r_y\leftarrow y$
​  for $j=1:n$
​    $r_a\leftarrow A(:,j)$
​    $r_y=r_y+r_a{r}_x(j)$
​  end
  $y\leftarrow r_y$

需要$(3+n)n$次读写操作，而对于外积计算：

​  $r_x\leftarrow x,r_y\leftarrow y$
​  for $j=1:n$
​    $r_a\leftarrow A(:,j)$
    $r_a=r_a+r_y{r}_x(j)$
    $A(:,j)\leftarrow r_a$
  end

需要$(2+2n)n$次读写操作，强弱悬殊高下立判。

3.3 存取优化-Gaxpy LU分解

Gaxpy的目标是结果以向量的形式来读写。那么这里计划的目标是在第$j$步时计算出$L$矩阵和$U$矩阵的第$j$列。当$j=1$时，由于$U(2:n,1)$的值皆为$0$，比较$A=LU$的第一列可知：
$$L(2:n,j)=A(2:n,1)/A(1,1),U(1,1)=A(1,1)$$
那么假设现在$L(:,1:j-1)$和$U(1:j-1,1:j-1)$已知（隐含$A$已知），那么根据以下算式可以求出$U(1:j-1,j)$：
$$A(1:j-1,j)=L(1:j-1,1:j-1)\cdot U(1:j-1,j)$$
当$U(1:j-1,j)$已求出之后，根据：
$$A(j:n,j)=L(j:n,1:j)\cdot U(1:j,j)$$
可以求出$U(j,j)$和$L(j+1:n,j)$。为了缩减Gaxpy LU分解所需要的空间，可以添加$n$维的辅助向量$v$:
$$v(j:n)=A(j:n,j)-L(j:n,1:j-1)\cdot U(1:j-1,j)$$
于是$L(j+1:n,j)=v(j+1:n)/v(j)$，$U(j,j)=v(j)$。算法伪代码如下：

  for $j=1:n$
    if $j=1$
      $v=A(:,1)$
    else
      $\tilde{a}=A(:,j)$
      Solve $L(1:j-1,1:j-1)\cdot z=\tilde{a}(1:j-1)$ for $z\in {\mathbb{R}}^{j-1}$
      $U(1:j-1,j)=z$
    end
    $U(j,j)=v(j)$
    $L(j+1:n,j)=v(j+1:n)/v(j)$
​  end

下面是C语言代码：

void gaxpyLU_decomposition(double* A, double* v, size_t n)
{
    size_t i, j, k;
    v[0] = A[0];
    for (i=1; i<n; i++)
    {
        v[i] = A[i * n];
        A[i * n] = A[i * n] / v[0];
    }
    for (i=1; i<n; i++)
    {
        for (j=1; j<i; j++)
            for (k=0; k<j; k++)
                A[j * n + i] -= A[j * n + k] * A[k * n + i];
        for (j=i; j<n; j++)
        {
            v[j] = A[j * n + i];
            for (k=0; k<i; k++)
                v[j] -= A[j * n + k] * A[k * n + i];
        }
        A[i * n + i] = v[i];
        for (j=i+1; j<n; j++)
            A[j * n + i] = v[j] / v[i];
    }
}

3.4 行方向Gaxpy LU分解（本文独有）

​ 根据内存中矩阵存储方式（按行或列），选择不同的分解方式可以减少数据存取的次数。下面推导按行方向的Gaxpy LU分解。

根据下角矩阵的性质可知：
$$L(i,j)=0,i,j\in {\mathbb{N}}^{+},i<j$$

$$U(i,j)=0,i,j\in {\mathbb{N}}^{+},i>j$$

计划目标是在第$j$步时计算出第$j$行的值（包括$L$和$U$），当$j=1$时可知：
$$L(1,1)=1,L(1,2:n)=0⟹U(1,i)=A(1,i),i=1,2,\dots ,n$$
现在假设$L(1:j-1,1:j-1)$和$U(1:j-1,1:j-1)$已知，那么通过
$$A(j,1:j-1)=L(j,1:j-1)\cdot U(1:j-1,1:j-1)$$
可以求出$L(j,1:j-1)，$已知$L(j,j)=1$，根据：
$$A(j,:)=L(j,1:j)\cdot U(1:j,:)$$
从而求出$U(j,:)$。该算法不需要辅助空间，可以进行原地分解，但是在伪代码中为了清晰还是区分了$L$和$U$矩阵，伪代码如下：

  $A(2:n,1)/=A(1,1)$
  for $j=2:n$
    Solve $z\cdot U(1:j-1,1:j-1)$ for $z\in {\mathbb{R}}^{j-1}$
    $L(j,1:j-1)=z,L(j,j)=1$
    Solve $L(j,1:j)\cdot w=A(j,:)$ for $w\in {\mathbb{R}}^j$
    $U(j,:)=w$
  end

C语言代码如下：

void rowbased_gaxpyLU_decomposition(double* A, size_t n)
{
    size_t i, j, k;
    double r;
    for (i=1; i<n; i++)
    {
    	A[i * n] /= A[0];
		for (j=1; j<i; j++)
		{
			for (k=0; k<j; k++) A[i * n + j] -= A[i * n + k] * A[k * n + j];
			A[i * n + j] /= A[j * n + j];
		}
		for (j=i; j<n; j++)
			for (k=0; k<i; k++) A[i * n + j] -= A[i * n + k] * A[k * n +j];
	}
}

4 数值稳定性和主元分解

4.1 稳定性分析

在$LU$分解后若出现数值极大的元素，会导致误差和数值不稳定，原因请参考《Matrix Computations》第4版的3.3节。下面举一个简单的例子解释：
$$A=\left[\begin{array}{cc}0.0001& 1\\ 1& 1\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 10000& 1\end{array}\right]\left[\begin{array}{cc}0.0001& 1\\ 0& -9999\end{array}\right]=LU\text{\hspace{1em}(4.1.1)}$$
$A$中数值极小的主元导致了分解后特别大数值的出现，克服此难题的一个方法是进行行交换，若使用行交换矩阵$P$：
$$P=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]$$
来左乘矩阵$A$之后再进行$LU$分解：
$$PA=\left[\begin{array}{cc}1& 1\\ 0.0001& 1\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 0.0001& 1\end{array}\right]\left[\begin{array}{cc}1& 1\\ 0& 0.9999\end{array}\right]=LU\text{\hspace{1em}(4.1.2)}$$
可以保证数值的稳定。

若要交换$A$矩阵的第$i$列和第$j$列，只要交换单位矩阵$I$的第$i$列和第$j$列得到矩阵$P$左乘$A$矩阵即可。

4.2 列主元$LU$分解

参考2.2节中的高斯消元法，为了保证数值的稳定，在第$j$步消元时，需要左乘行交换矩阵${\Pi }_i$将第$j$列中从$j$到$n$位置值最大的元素交换到第$j$行，再进行高斯消元，完整的步骤如下：
$$M_{n-1}{\Pi }_{n-1}\cdots M_1{\Pi }_{1}A=U\text{\hspace{1em}(4.2.1)}$$
接下来需要找出$L$矩阵，从最简单的关系开始：
$$PA=LU,P={\Pi }_{n-1}\cdots {\Pi }_1\text{\hspace{1em}(4.2.2)}$$
结合式（4.2.1）可以得出以下结果：
$${\tilde{M}}_{n-1}{\tilde{M}}_{n-2}\cdots {\tilde{M}}_{1}PA=U\text{\hspace{1em}(4.2.3)}$$
由于${\Pi }_{i}A{\Pi }_i=A$，显然：
$${\tilde{M}}_k=({\Pi }_{n-1}\cdots {\Pi }_{k+1})M_k({\Pi }_{k+1}\cdots {\Pi }_{n-1})\text{\hspace{1em}(4.2.4)}$$
由线性代数中所学的知识可知，左乘矩阵是进行行变换，而右乘矩阵则是进行列变换。直观地说，交换单位矩阵$I$的第$i$行和第$j$行得到行交换矩阵$\Pi$，用$\Pi$左乘任意矩阵的结果便是交换该矩阵的第$i$行和第$j$行。
在进行$LU$分解时为保持数值稳定，会在每一步的高斯消元时进行行交换，并保存行交换的信息在矩阵$P$中。最终分解的公式变为：
$$PA=LU\text{\hspace{1em}(4.2.5)}$$
以下算法伪代码完成了$4.2.5$式的计算：
   $P=I$
   for $k=1:n-1$
     Find $\mu$ where $|A(\mu ,k)|=\Vert A(k:n,k){\Vert }_{\infty }$
     $Exchange(A(\mu ,:),A(k,:))$
     $Exchange(P(\mu ,:),P(k,:))$
     if $A(k,k)\ne 0$
       $\rho =k+1:n$
       $A(\rho ,k)=A(\rho ,k)/A(k,k)$
       $A(\rho ,\rho )=A(\rho ,\rho )-A(\rho ,k)A(k,\rho )$
     end
   end