Main Content

lu

LU 矩阵分解

全页折叠

语法

[L,U] = lu(A)
[L,U,P] = lu(A)
[L,U,P] = lu(A,outputForm)
[L,U,P,Q] = lu(S)
[L,U,P,Q,D] = lu(S)
[___] = lu(S,thresh)
[___] = lu(___,outputForm)

说明

示例

[L,U] = lu(A) 将满矩阵或稀疏矩阵 A 分解为一个上三角矩阵 U 和一个经过置换的下三角矩阵 L，使得 A = L*U。

示例

[L,U,P] = lu(A) 还返回一个置换矩阵 P，并满足 A = P'*L*U。在此语法中，L 是单位下三角矩阵，U 是上三角矩阵。

示例

[L,U,P] = lu(A,outputForm) 以 outputForm 指定的格式返回 P。将 outputForm 指定为 'vector' 会将 P 返回为一个置换向量，并满足 A(P,:) = L*U。

示例

[L,U,P,Q] = lu(S) 将稀疏矩阵 S 分解为一个单位下三角矩阵 L、一个上三角矩阵 U、一个行置换矩阵 P 以及一个列置换矩阵 Q，并满足 P*S*Q = L*U。

[L,U,P,Q,D] = lu(S) 还返回一个对角缩放矩阵 D，并满足 P*(D\S)*Q = L*U。行缩放通常会使分解更为稀疏和稳定。

[___] = lu(S,thresh) 可结合上述任意输出参数组合指定 lu 使用的主元消去策略的阈值。根据指定的输出参数的数量，对 thresh 输入的要求及其默认值会有所不同。有关详细信息，请参阅 thresh 参数描述。

示例

[___] = lu(___,outputForm) 以 outputForm 指定的格式返回 P 和 Q。将 outputForm 指定为 'vector' 以将 P 和 Q 返回为置换向量。您可以使用上述语法中的任何输入参数组合。

示例

全部折叠

打开实时脚本

计算矩阵的 LU 分解并检查生成的因子。通过 LU 分解，可以将一个矩阵 A 分解为一个上三角矩阵 U、一个下三角矩阵 L 和一个置换矩阵 P，并满足 PA=LU。这些矩阵描述了对矩阵执行高斯消去法，直至其化简成为行阶梯形矩阵所需的步骤。L 矩阵包含所有乘数，置换矩阵 P 解释行交换。

创建一个 3×3 矩阵并计算 LU 因子。

A = [10 -7 0
     -3  2 6
      5 -1 5];
[L,U] = lu(A)
L = 3×3

    1.0000         0         0
   -0.3000   -0.0400    1.0000
    0.5000    1.0000         0


U = 3×3

   10.0000   -7.0000         0
         0    2.5000    5.0000
         0         0    6.2000

将因子相乘以重新创建 A。在双输入语法的基础上，lu 将置换矩阵 P 直接合并到 L 因子中，使得返回的 L 实际上等于 P'*L，因此 A = L*U。 

L*U
ans = 3×3

    10    -7     0
    -3     2     6
     5    -1     5

您可以指定三个输出以将置换矩阵与 L 中的乘数分开。

[L,U,P] = lu(A)
L = 3×3

    1.0000         0         0
    0.5000    1.0000         0
   -0.3000   -0.0400    1.0000


U = 3×3

   10.0000   -7.0000         0
         0    2.5000    5.0000
         0         0    6.2000


P = 3×3

     1     0     0
     0     0     1
     0     1     0


P'*L*U
ans = 3×3

    10    -7     0
    -3     2     6
     5    -1     5
打开实时脚本

通过执行 LU 分解并使用因子来简化问题，对线性系统求解。使用反斜杠运算符和 decomposition 对象将结果与其他方法进行比较。

创建一个 5×5 幻方矩阵并求解线性系统 Ax=b，其中 b 的所有元素等于 65，即幻数和。由于 65 是此矩阵的幻数和（各行之和与各列之和均为 65），因此 x 的预期解是由 1 组成的向量。

A = magic(5);
b = 65*ones(5,1);
x = A\b
x = 5×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000

对于泛型方阵，反斜杠运算符使用 LU 分解计算线性系统的解。LU 分解将 A 表示为三角矩阵的乘积，它可以通过代换公式轻松求解涉及三角矩阵的线性系统。

要重新创建由反斜杠运算符计算的答案，请计算 A 的 LU 分解。然后，使用因子来求解两个三角线性系统：

y = L\(P*b);
x = U\y;

这种在求解线性系统之前预先计算矩阵因子的方法可以在求解许多线性系统时提高性能，因为分解仅发生一次而不需要重复。

[L,U,P] = lu(A)
L = 5×5

    1.0000         0         0         0         0
    0.7391    1.0000         0         0         0
    0.4783    0.7687    1.0000         0         0
    0.1739    0.2527    0.5164    1.0000         0
    0.4348    0.4839    0.7231    0.9231    1.0000


U = 5×5

   23.0000    5.0000    7.0000   14.0000   16.0000
         0   20.3043   -4.1739   -2.3478    3.1739
         0         0   24.8608   -2.8908   -1.0921
         0         0         0   19.6512   18.9793
         0         0         0         0  -22.2222


P = 5×5

     0     1     0     0     0
     1     0     0     0     0
     0     0     0     0     1
     0     0     1     0     0
     0     0     0     1     0


y = L\(P*b);
x = U\y
x = 5×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000

在使用专用分解来求解线性系统时 decomposition 对象也很有用，因为预先计算矩阵因子可获得许多性能优势，而不需要知道如何使用这些因子。使用 'lu' 类型的分解对象重新创建相同的结果。

dA = decomposition(A,'lu');
x = dA\b
x = 5×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
打开实时脚本

计算稀疏矩阵的 LU 分解并验证恒等式 L*U = P*S*Q。

基于 Buckminster-Fuller 多面穹顶连接图创建一个 60×60 稀疏邻接矩阵。

S = bucky;

使用带有四个输出的稀疏矩阵语法计算 S 的 LU 分解，以返回行和列置换矩阵。

[L,U,P,Q] = lu(S);

使用 P*S*Q 置换 S 的行和列，并将结果与三角因子 L*U 的乘积进行比较。二者之差的 1-范数在舍入误差之内，这表明 L*U = P*S*Q。

e = P*S*Q - L*U;
norm(e,1)
ans = 2.2204e-16
打开实时脚本

计算矩阵的 LU 分解。通过将行置换返回为向量（而不是矩阵）来节省内存。

创建一个 1000×1000 随机矩阵。

A = rand(1000);

计算 LU 分解并将置换信息存储为矩阵 P。将其结果与将置换信息存储为向量 p 的分解结果进行比较。矩阵越大，使用置换向量的内存效率越高。

[L1,U1,P] = lu(A);
[L2,U2,p] = lu(A,'vector');
whos P p
  Name         Size                Bytes  Class     Attributes

  P         1000x1000            8000000  double              
  p            1x1000               8000  double

使用置换向量还可以节省后续操作中的执行时间。例如，您可以使用前面介绍的 LU 分解来求解线性系统 Ax=b。尽管从置换向量和置换矩阵求得的解是相等的（取决于舍入），但使用置换向量求解所需的时间通常略少。

打开实时脚本

将在使用和不使用列置换的情况下计算稀疏矩阵的 LU 分解的结果进行比较。

加载 west0479 矩阵，这是一个实数值 479×479 稀疏矩阵。

load west0479
A = west0479;

通过调用带三个输出的 lu 来计算 A 的 LU 分解。生成 L 因子和 U 因子的 spy 图。

[L,U,P] = lu(A);
subplot(1,2,1)
spy(L)
title('L factor')
subplot(1,2,2)
spy(U)
title('U factor')

Figure contains 2 axes objects. axes object 1 with title L factor, xlabel nz = 10351 contains a line object which displays its values using only markers. axes object 2 with title U factor, xlabel nz = 6046 contains a line object which displays its values using only markers.

现在，使用带四个输出的 lu 计算 A 的 LU 分解，这将置换 A 的列以减少因子中的非零元素数。如果不使用列置换，则产生的因子更为稀疏。

[L,U,P,Q] = lu(A);
subplot(1,2,1)
spy(L)
title('L factor')
subplot(1,2,2)
spy(U)
title('U factor')

Figure contains 2 axes objects. axes object 1 with title L factor, xlabel nz = 1854 contains a line object which displays its values using only markers. axes object 2 with title U factor, xlabel nz = 2391 contains a line object which displays its values using only markers.

输入参数

全部折叠

输入矩阵。A 可以是满矩阵或稀疏矩阵，也可以是方阵或矩形矩阵。

数据类型： single | double
复数支持： 是

稀疏输入矩阵。S 可以是方阵或矩形矩阵。

数据类型： double
复数支持： 是

稀疏矩阵的主元阈值，指定为标量或二元素向量。有效值介于区间 [0 1] 之间。指定 thresh 的方式取决于在调用 lu 时所指定的输出数量：

  • 当输出为三个或少于三个时，thresh 必须是标量，默认值为 1.0。

  • 当输出为四个或多于四个时，thresh 可以是标量或二元素向量。默认值为 [0.1 0.001]。如果将 thresh 指定为标量，则仅替换向量中的第一个值。

从较高层面讲，此输入使您能够在精确度和总执行时间之间进行权衡。thresh 的值越小，LU 因子越稀疏，但解可能不准确。值越大，生成的解越准确（并非总是如此），并且总工作量和内存使用率越高。

lu 首先根据输出参数的数量、然后根据待分解矩阵的属性选择主元消去策略。在任何情况下，将阈值设置为 1.0 都会导致部分主元消去，而将它们设置为 0 则会导致仅根据生成的矩阵的稀疏性选择主元。L 的所有值的绝对值都小于或等于 1/min(thresh)。

  • 三个或少于三个的输出参数 - 如果算法满足以下等式，算法会选择对角主元

    A(j,j) >= thresh * max(abs(A(j:m,j)))
    否则，它会选择包含绝对值最大的元素所在的行。

  • 对称主元消去策略 - 如果 S 是大体结构对称和对角线上元素多为非零元素的稀疏方阵，则 lu 将使用对称主元消去策略。对于此策略，如果算法满足以下不等式，则算法会选择对角主元 j：

    A(i,j) >= thresh(2) * max(abs(A(j:m,j)))
    如果对角线元素未通过此测试，则 lu 会选择满足以下不等式的最稀疏的行 i：
    A(i,j) >= thresh(1) * max(abs(A(j:m,j)))

  • 非对称主元消去策略 - 如果 S 不满足对称主元消去策略的要求，则 lu 将使用非对称策略。在这种情况下，lu 会选择满足以下不等式的最稀疏的行 i：

    A(i,j) >= thresh(1) * max(abs(A(j:m,j)))
    对于 thresh(1)，值为 1.0 会执行传统的部分主元消去。L 中的元的绝对值都小于或等于 1/thresh(1)。当使用非对称策略时，不会使用 thresh 输入向量的第二个元素。

注意

在极少数情况下，错误的分解会导致 P*S*Q ≠ L*U。如果发生这种情况，请将 thresh 增至最大值 1.0（常规部分主元消去），并重试。

置换输出的形状，指定为 'matrix' 或 'vector'。此标志控制 lu 是将行置换结果 P 和列置换结果 Q 返回为置换矩阵还是返回为置换向量。

如果返回为矩阵，则输出 P 和 Q 满足以下恒等式：

  • 三个输出 - P 满足 P*A = L*U。

  • 四个输出 - P 和 Q 满足 P*S*Q = L*U。

  • 五个输出 - P、Q 和 D 满足 P*(D\S)*Q = L*U。

如果返回为向量，则输出 P 和 Q 满足以下恒等式：

  • 三个输出 - P 满足 A(P,:) = L*U

  • 四个输出 - P 和 Q 满足 S(P,Q) = L*U

  • 五个输出 - P、Q 和 D 满足 D(:,P)\S(:,Q) = L*U。

示例： [L,U,P] = lu(A,'vector')

输出参数

全部折叠

下三角因子，以矩阵形式返回。L 的形式取决于是否在一个单独的输出中返回行置换结果 P：

  • 如果指定第三个输出 P，则 L 返回为一个单位下三角矩阵（即，主对角线上元素为 1 的下三角矩阵）。

  • 如果未指定第三个输出 P，则 L 返回为一个单位下三角矩阵的行置换矩阵。具体而言，它是三输出情况下所返回输出 P 和 L 的 P'*L 的乘积。

上三角因子，以上三角矩阵形式返回。

行置换，以置换矩阵形式返回；如果指定 'vector' 选项，则以置换向量形式返回。使用此输出可提高计算的数值稳定性。

有关此输出满足的恒等式的说明，请参阅 outputForm。

列置换，以置换矩阵形式返回；如果指定 'vector' 选项，则以置换向量形式返回。使用此输出可减少稀疏矩阵因子矩阵中的填充量（非零元素数）。

有关此输出满足的恒等式的说明，请参阅 outputForm。

行缩放，以对角矩阵形式返回。D 用于缩放 S 中的值，使得 P*(D\S)*Q = L*U。行缩放通常（但不是始终）会使分解更为稀疏和稳定。

算法

使用高斯消去法的变体计算 LU 分解。计算精确解取决于原始矩阵 cond(A) 的条件数的值。如果矩阵具有较大的条件数（接近奇异矩阵），则计算的分解可能不准确。

LU 分解是使用 inv 得到逆矩阵和使用 det 得到行列式的关键步骤。同时，它还是使用运算符 \ 和 / 求线性方程解或矩阵除法的基础。这意味着 lu 的这些数值限制也会存在于依赖它的这些函数中。

参考

[1] Gilbert, John R., and Tim Peierls. “Sparse Partial Pivoting in Time Proportional to Arithmetic Operations.” SIAM Journal on Scientific and Statistical Computing 9, no. 5 (September 1988): 862–874. https://doi.org/10.1137/0909058.

[2] Anderson, E., ed. LAPACK Users’ Guide. 3rd ed. Software, Environments, Tools. Philadelphia: Society for Industrial and Applied Mathematics, 1999. https://doi.org/10.1137/1.9780898719604.

[3] Davis, Timothy A. "Algorithm 832: UMFPACK V4.3 – an unsymmetric-pattern multifrontal method." ACM Transactions on Mathematical Software 30, no. 2 (June 2004): 196–199. https://doi.org/10.1145/992200.992206.

扩展功能

版本历史记录

在 R2006a 之前推出

另请参阅

| | | | | | |

主题