Matlab快速入门之线性代数:矩阵分解

4.3
(3)

本文属于Matlab快速入门之线性代数的第三篇,即矩阵分解,主要包括:Cholesky 分解、LU 分解、QR 分解、对分解使用多线程计算。讨论的这三种矩阵分解利用了三角形矩阵,其中对角线上下的所有元素都为零。涉及三角矩阵的线性方程组可以使用前代或回代方法轻松快捷地求解。

Matlab快速入门之线性代数:矩阵分解

Cholesky 分解

Cholesky 分解将对称矩阵表示为三角矩阵与其转置的积A = R′R,其中,R 是上三角矩阵。

并非所有对称矩阵都可以通过这种方式进行分解;采用此类分解的矩阵被视为正定矩阵。这表明,A 的所有对角线元素都是正数,并且非对角线元素“不太大”。帕斯卡矩阵提供了有趣的示例。在本章中,示例矩阵 A 为 3×3 帕斯卡矩阵,暂时转换为 6×6,A 的元素为二项式系数。每个元素都是其北向和西向邻点之和。Cholesky 分解可以使用chol函数获得:

>> A = pascal(6)

A =

     1     1     1     1     1     1
     1     2     3     4     5     6
     1     3     6    10    15    21
     1     4    10    20    35    56
     1     5    15    35    70   126
     1     6    21    56   126   252

>> R = chol(A)

R =

     1     1     1     1     1     1
     0     1     2     3     4     5
     0     0     1     3     6    10
     0     0     0     1     4    10
     0     0     0     0     1     5
     0     0     0     0     0     1
Matlab快速入门之线性代数:矩阵分解

这些元素同样为二项式系数。R'*R 等于 A 的情况说明了涉及二项式系数的积之和的单位矩阵。

注意

Cholesky 分解也适用于复矩阵。采用 Cholesky 分解的任何复矩阵都满足
A′ = A,并且被视为 Hermitian 正定矩阵

通过 Cholesky 分解,可以将线性方程组Ax = b替换为R′Rx = b

由于反斜杠运算符能识别三角形方程组,因此这可以在 MATLAB 环境中通过以下表达式快速进行求解:

x = R\(R'\b)

如果 A 为 n×n,则 chol(A) 的计算复杂度为 O(n3),但后续的反斜杠解的复杂度仅为 O(n2)。

LU 分解

LU 分解(或高斯消去法)将任何方阵 A 都表示为下三角矩阵和上三角矩阵的置换之积:A = LU

其中,L 是对角线元素为 1 的下三角矩阵的置换,U 是上三角矩阵。

出于理论和计算原因,必须进行置换。矩阵

[0  1
1  0]

在不交换其两行的情况下不能表示为三角矩阵的积。尽管矩阵

[ε  1
1   0]

可以表示为三角矩阵之积,但当ε很小时,因子中的元素也会很大并且会增大误差,因此即使置换并非完全必要,它们也是所希望的。部分主元消去法可确保 L 的元素的模以 1 为限,并且 U 的元素并不大于 A 的元素。

例如:

[L,U] = lu(B)

L =
    1.0000         0         0
    0.3750    0.5441    1.0000
    0.5000    1.0000         0

U =
    8.0000    1.0000    6.0000
         0    8.5000   -1.0000
         0         0    5.2941

通过对 A 执行 LU 分解,可以

A*x = b

使用以下表达式快速对线性方程组求解

x = U\(L\b)

行列式和逆矩阵是通过 LU 分解使用以下表达式进行计算的

det(A) = det(L)*det(U)

以及

inv(A) = inv(U)*inv(L)

也可以使用 det(A) = prod(diag(U)) 计算行列式,但行列式的符号可能会相反。

QR 分解

正交矩阵或包含正交列的矩阵为实矩阵,其列全部具有单位长度并且相互垂直。如果 Q 为正交矩阵,则

QTQ = I

其中 I 是单位矩阵。

最简单的正交矩阵为二维坐标旋转:

[ cos(θ)  sin(θ)
 -sin(θ)  cos(θ)]

对于复矩阵,对应的术语为单位。由于正交矩阵和单位矩阵会保留长度、保留角度并且不会增大误差,因此适用于数值计算。

正交或 QR 分解将任何矩形矩阵表示为正交或酉矩阵与上三角矩阵的积。此外,也可能涉及列置换。

A = QR

AP = QR

其中,Q 为正交或单位矩阵,R 为上三角矩阵,P 为置换向量。

QR 分解有四种变化形式 – 完全大小或合适大小,以及使用列置换或不使用列置换。

超定线性方程组涉及行数超过列数的矩形矩阵,也即 m×n 并且 m > n。完全大小的 QR 分解可生成一个方阵(m×m 正交矩阵 Q)和一个矩形 m×n 上三角矩阵 R:

>> C=gallery('uniformdata',[5 4], 0);
>> [Q,R] = qr(C)

Q =

    0.6191    0.1406   -0.1899   -0.5058    0.5522
    0.1506    0.4084    0.5034    0.5974    0.4475
    0.3954   -0.5564    0.6869   -0.1478   -0.2008
    0.3167    0.6676    0.1351   -0.1729   -0.6370
    0.5808   -0.2410   -0.4695    0.5792   -0.2207


R =

    1.5346    1.0663    1.2010    1.4036
         0    0.7245    0.3474   -0.0126
         0         0    0.9320    0.6596
         0         0         0    0.6648
         0         0         0         0
Matlab快速入门之线性代数:矩阵分解

在许多情况下,Q 的最后 m – n 列是不需要的,因为这些列会与 R 底部的零相乘。因此,精简 QR 分解可生成一个矩形矩阵(具有正交列的 m×n Q)以及一个方阵 n×n 上三角矩阵 R。对于 5×4 示例,不会节省太多内存,但是对于更大的大量矩形矩阵,在时间和内存方面的节省可能会很重要:

>> [Q,R] = qr(C,0)	
Q =

    0.6191    0.1406   -0.1899   -0.5058
    0.1506    0.4084    0.5034    0.5974
    0.3954   -0.5564    0.6869   -0.1478
    0.3167    0.6676    0.1351   -0.1729
    0.5808   -0.2410   -0.4695    0.5792


R =

    1.5346    1.0663    1.2010    1.4036
         0    0.7245    0.3474   -0.0126
         0         0    0.9320    0.6596
         0         0         0    0.6648
Matlab快速入门之线性代数:矩阵分解

与 LU 分解相比,QR 分解不需要进行任何消元或置换。但是,可选的列置换(因存在第三个输出参数而触发)对检测奇异性或秩亏是很有帮助的。在分解的每一步,未分解的剩余矩阵的列(范数最大)将用作该步骤的基础。这可以确保 R 的对角线元素以降序排列,并且各列之间的任何线性相关性肯定能够通过检查这些元素来显示。对于此处提供的小示例,C 的第二列的范数大于第一列的范数,因此这两列被交换:

>> [Q,R,P] = qr(C)

Q =
   -0.3522    0.8398   -0.4131
   -0.7044   -0.5285   -0.4739
   -0.6163    0.1241    0.7777

R =
  -11.3578   -8.2762
         0    7.2460
         0         0

P =
     0     1
     1     0

组合了合适大小和列置换后,第三个输出参数为置换向量而不是置换矩阵:

>> [Q,R,p] = qr(C,0)

Q =
   -0.3522    0.8398
   -0.7044   -0.5285
   -0.6163    0.1241

R =
  -11.3578   -8.2762
         0    7.2460


p =
     2     1

QR 分解可将超定线性方程组转换为等效的三角形方程组。表达式

norm(A*x - b)

等于

norm(Q*R*x - b)

与正交矩阵相乘可保留欧几里德范数,因此该表达式也等于

norm(R*x - y)

其中 y = Q'*b。由于 R 的最后 m-n 行为零,因此该表达式将分为两部分:

norm(R(1:n,1:n)*x - y(1:n))

以及

norm(y(n+1:m))

如果 A 具有满秩,则可以对 x 求解,使这些表达式中的第一个表达式为零。然后,第二个表达式便可以提供残差范数。如果 A 没有满秩,则可以通过 R 的三角形结构体对最小二乘问题求基本解。

对分解使用多线程计算

对于许多线性代数函数和按元素的数值函数,MATLAB 软件支持多线程计算。这些函数将自动在多个线程上执行。要使函数或表达式在多个 CPU 上更快地执行,必须满足许多条件:

  1. 函数执行的运算可轻松划分为并发执行的多个部分。这些部分必须能够在进程之间几乎不通信的情况下执行。它们应需要很少的序列运算。
  2. 数据大小足以使并发执行的任何优势在重要性方面超过对数据分区和管理各个执行线程所需的时间。例如,仅当数组包含数千个或以上的元素时,大多数函数才会加速。
  3. 运算未与内存绑定;处理时间不受内存访问时间控制。一般而言,复杂函数比简单函数速度更快。

对于大型双精度数组(约 10,000 个元素),lu 和 qr 会大幅增加速度。

共计3人评分,平均4.3

到目前为止还没有投票~

很抱歉,这篇文章对您没有用!

让我们改善这篇文章!

告诉我们我们如何改善这篇文章?

文章目录

转载文章,原文出处:MathWorks官网,由古哥整理发布

如若转载,请注明出处:https://iymark.com/articles/3260.html

(0)
微信公众号
古哥的头像古哥管理团队
上一篇 2022年09月23日 20:08
下一篇 2022年09月25日 17:49

你可能感兴趣的文章

发表回复

登录后才能评论
微信小程序
微信公众号