LU分解(C++)

LU分解是一种重要的数值线性代数技术, 用于解决线性方程组和矩阵求逆等问题. 在科学工程领域, 经常需要解决形如$Ax=b$的线性方程组, 其中$A$是系数矩阵, $x$是未知向量, $b$是已知向量. LU分解是一种将系数矩阵$A$分解为一个下三角矩阵$L$和一个上三角矩阵$U$的方法, 即$A=LU$. 这个分解有许多优点, 其中之一是它可以帮助我们更有效地解决多个不同的右端向量$b$对应的线性方程组, 而无需每次都重新分解$A$. 此外, LU分解也有助于计算矩阵的逆, 因为一旦我们得到了$A=LU$的分解, 就可以相对容易地计算出$A$的逆. 因此, LU分解是许多数值计算和线性代数问题的基础.

LU分解在数值计算中具有广泛的应用, 包括但不限于以下几个方面:

1. 线性方程组求解 LU分解可用于有效解决多个不同右端向量的线性方程组, 减少了重复分解系数矩阵的计算开销.
2. 矩阵求逆 一旦得到$A=LU$的分解, 可以相对容易地计算矩阵$A$的逆矩阵. 这在各种科学计算和工程应用中非常有用.
3. 数值稳定性 LU分解可以帮助分析和改善数值算法的稳定性, 以减小舍入误差的影响.
4. 最小二乘法 LU分解可用于最小二乘拟合等问题, 其中需要求解超定或欠定线性方程组.

本文仅考虑不选主元素的三角分解方法, 即$A$的顺序主子式都不等于0, 否则需要对矩阵$A$左乘一个排列矩阵$P$.

算法描述

LU分解的数学原理基于Gauss消元法. 它的核心思想是将系数矩阵$A$通过一系列行变换变成一个上三角矩阵$U$, 同时记录下每次行变换的过程, 以构造下三角矩阵$L$.

设$A$是一个$n\times n$的矩阵, LU分解的目标是找到下三角矩阵$L$和上三角矩阵$U$, 使得$A=LU$. 具体步骤包括：

1. 将$L$初始化为单位下三角矩阵, 将$U$初始化为$A$的副本.
2. 针对第一列, 使用行变换操作将$U$的第一列元素变成零, 同时记录行变换操作到$L$的第一列.
3. 重复上述步骤, 依次处理第二列, 第三列, 直到处理完最后一列, 得到完整的LU分解.

数学上, 行变换操作是通过矩阵乘法来表示的, 这些操作将$A$的行变换为$U$的行, 同时更新$L$. 最终, 得到的矩阵$U$就是原矩阵$A$的上三角部分, 而矩阵$L$则包含了所有行变换的信息.

实际上, 如果每次都进行消元, 可能导致不稳定的情形出现, 且效率较低. 我们可以直接从LU分解的结论$A=LU$出发, 利用矩阵乘法的定义, 我们可以得到$n^2$个方程, 通过化简计算, 可得如下快速计算LU分解的算法:

设$A=(a_{ij}),L=(l_{ij}),U=(u_{ij})$, 首先由$A$的第1行第1列可以计算出$U$的第1行和$L$的第1列:

$$u_{1j}=a_{1j},j=1,2,?,n$$

$$l_{k1}=\frac{a_{k1}}{u_{11}},k=2,3,?,n$$

下面假设$U$的$1$到$k?1$行, $L$的$1$到$k?1$列均已算出, 则有:

$$u_{kj}=a_{kj}?\sum _{r=1}^{k?1}{l}_{kr}{u}_{rj},j=k,k+1,?,n$$

$$l_{ik}=\frac{1}{u_{kk}}\left(a_{ik}?\sum _{r=1}^{k?1}{l}_{ir}{u}_{rk}\right),i=k,k+1,?,n$$

根据如上递推公式即可算出$L$和$U$的全部元素.

算法实现

#include <armadillo>
using namespace arma;
/*
 * LU分解
 * L:下三角矩阵
 * U:上三角矩阵
 * A:待分解矩阵
 * e:精度
 *
 * 返回(bool):
 *  true : 分解失败
 *  false: 分解成功
 */
bool LU(mat &L, mat &U, const mat &A, const double &e = 1e-6)
{
    if (A.n_cols == 1)
    {
        L.ones(1, 1);
        U.resize(1, 1);
        if (abs(U.at(0, 0) = A.at(0, 0)) < e)
            return true;
        return false;
    }
    L.eye(A.n_cols, A.n_cols);
    U.zeros(A.n_cols, A.n_cols);
    unsigned n(A.n_cols - 1);
    for (unsigned i(0); i != n; ++i)
    {
        U.at(i, i) = A.at(i, i);
        for (unsigned k(0); k != i; ++k)
            U.at(i, i) -= L.at(i, k) * U.at(k, i);
        if (abs(U.at(i, i)) < e)
            return true;
        for (unsigned j(i + 1); j != A.n_cols; ++j)
        {
            L.at(j, i) = A.at(j, i);
            U.at(i, j) = A.at(i, j);
            for (unsigned k(0); k != i; ++k)
            {
                U.at(i, j) -= L.at(i, k) * U.at(k, j);
                L.at(j, i) -= L.at(j, k) * U.at(k, i);
            }
            L.at(j, i) /= U.at(i, i);
        }
    }
    U.at(n, n) = A.at(n, n);
    for (unsigned i(0); i != n; ++i)
        U.at(n, n) -= L.at(n, i) * U.at(i, n);
    if (abs(U.at(n, n)) < e)
        return true;
    return false;
}

实际上armadillo库提供了lu函数, 也可以直接使用arma::lu.

实例分析

对以下矩阵进行LU分解:

$$A=\left(\begin{array}{cccc}6& 2& 1& ?1\\ 2& 4& 1& 0\\ 1& 1& 4& ?1\\ ?1& 0& ?1& 3\end{array}\right)$$

代入程序求得

$$L=\left(\begin{array}{cccc}1.0000& 0& 0& 0\\ 0.3333& 1.0000& 0& 0\\ 0.1667& 0.2000& 1.0000& 0\\ ?0.1667& 0.1000& ?0.2432& 1.0000\end{array}\right)$$

$$U=\left(\begin{array}{cccc}6.0000& 2.0000& 1.0000& ?1.0000\\ 0& 3.3333& 0.6667& 0.3333\\ 0& 0& 3.7000& ?0.9000\\ 0& 0& 0& 2.5811\end{array}\right)$$