LU分解(C++)

原创于 2024-01-16 00:01:41 发布 · 2k 阅读

19 ·

CC 4.0 BY-SA版权

文章标签：

#c++ #算法

基于armadillo库的数值分析实践专栏收录该内容

17 篇文章

订阅专栏

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

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

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

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

算法描述

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

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

将 $L$ 初始化为单位下三角矩阵, 将 $U$ 初始化为 $A$ 的副本.
针对第一列, 使用行变换操作将 $U$ 的第一列元素变成零, 同时记录行变换操作到 $L$ 的第一列.
重复上述步骤, 依次处理第二列, 第三列, 直到处理完最后一列, 得到完整的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,\cdots,n$

$l_{k1}=\frac{a_{k1}}{u_{11}},k=2,3,\cdots,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,\cdots,n$

$l_{ik}=\frac1{u_{kk}}\left(a_{ik}-\sum_{r=1}^{k-1}l_{ir}u_{rk}\right),i=k,k+1,\cdots,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;
}