跳至主要内容

机器学习中的 QR 分解:详细指南

了解 QR 分解,这是一种将矩阵 A 分解为正交矩阵 Q 与上三角矩阵 R 之积的矩阵分解技术。
· 12分钟 阅读

使用 AI 探索

QR 分解是一种基础的矩阵分解技术,广泛应用于数据科学和机器学习的诸多领域。它能将复杂矩阵拆解为更简单的部分,因此在求解线性代数问题与提升算法稳定性方面尤为重要。

无论您是在进行普通最小二乘线性回归,还是用 Python 或 PySpark构建推荐引擎,理解 QR 分解都能显著增强您处理大型数据集与复杂数学模型的能力。

什么是 QR 分解?

对矩阵进行因式分解(分解)是指将该矩阵表示为两个或多个矩阵的乘积。换言之,矩阵分解允许您把一个矩阵拆解为更简单的组件,这些组件在进行矩阵乘法后可还原为原矩阵。我们之所以重视矩阵分解,是因为它能简化许多线性代数问题,使之更易求解。

QR 分解是一种特定的矩阵分解方法,用于将矩阵分解为正交矩阵Q与上三角矩阵R的乘积。我们常见的等式为 A = QR。

QR 分解可视化

QR 分解可视化。来源:ResearchGate

要理解 QR 分解,您需要注意在线性代数中,矩阵相乘并不遵循标量乘法的直觉规则。具体来说,矩阵乘法不满足交换律,一般 AB ≠ BA,因此“逆向”的概念也更为复杂。

Q 与 R 矩阵是什么?

我们来看看两个关键的矩阵:Q 和 R。

正交矩阵 Q

正交矩阵 Q 的特征是其转置也是其逆矩阵,数学上表示为 QTQ = QQT = I,其中 I 为单位矩阵。作个背景说明,转置是指沿着对角线翻转矩阵。

正交矩阵还具有其列向量构成一组正交归一基的性质。这意味着两点:其一,列向量两两正交,即矩阵中任意两个不同列的点积为零;其二,每一列的模(单位长度)为一。

正交矩阵之所以易于使用,是因为它能保持向量的长度及其点积不变。因此,如果您将一个向量与正交矩阵相乘,该向量的长度保持不变。同样地,如果您将两个向量同时与同一个正交矩阵相乘,它们之间的点积也保持不变。最后,由于正交矩阵保持长度和点积,变换后不会扭曲向量的几何性质。也就是说,如果两个向量在变换前是互相垂直的,那么在应用正交矩阵后它们仍将互相垂直。

上三角矩阵 R

上三角矩阵 R 的定义是主对角线以下的所有元素都为零。这种形式可以简化矩阵方程,特别适用于通过回代或其他方法求解线性方程组。在 QR 分解中,我们也可以说 R 表征了原矩阵的大小与方向信息。

执行 QR 分解的方法

我们来快速浏览几种进行 QR 分解的方法。

Gram-Schmidt 过程

Gram-Schmidt 过程是一种将一组线性无关向量系统地转换为正交归一集合的算法。该过程专门用于构造 QR 分解中的 Q 矩阵,在得到 Q 之后,我们可以推导出 R。

Givens 旋转

Givens 旋转是帮助进行 QR 分解的另一种方法。通过一系列平面旋转,将矩阵 A 的对角线以下元素消为零。每次旋转作用于一对坐标,既保持向量长度,又将矩阵变换为上三角形式。该方法对稀疏矩阵效率较高,因为每一步仅修改矩阵的一小部分。

Householder 反射

Householder 反射是一种更通用的技术,使用称为 Householder 矩阵的反射矩阵,系统性地消去矩阵 A 的次对角线以下元素。每个 Householder 矩阵将向量关于与所选向量正交的平面进行反射,从而把某列对角线以下元素归零。结果得到正交矩阵 Q 和上三角矩阵 R。Householder 反射因其数值稳定性而尤其适用于稠密矩阵。

在 R 中实现 QR 分解

本节中,我们将通过使用内置的 qr() 函数快速实现一次 QR 分解,然后再从零开始用 Gram-Schmidt 过程走一遍,以加深直觉。首先创建 matrix_A,两部分都会用到它。

matrix_A <- matrix(c(21, 21, 22.8, 21.4, 18.7, 18.1, 14.3, 24.4, 22.8, 19.2, 17.8, 16.4, 17.3, 15.2, 10.4, 10.4, 14.7, 32.4, 30.4, 33.9, 21.5, 15.5, 15.2, 13.3, 19.2, 27.3, 26, 30.4, 15.8, 19.7, 15, 21.4, 
              2.62, 2.875, 2.32, 3.215, 3.44, 3.46, 3.57, 3.19, 3.15, 3.44, 3.44, 4.07, 3.73, 3.78, 5.25, 5.424, 5.345, 2.2, 1.615, 1.835, 2.465, 3.52, 3.435, 3.84, 3.845, 1.935, 2.14, 1.513, 3.17, 2.77, 3.57, 2.78, 
              110, 110, 93, 110, 175, 105, 245, 62, 95, 123, 123, 180, 180, 180, 205, 215, 230, 66, 52, 65, 97, 150, 150, 245, 175, 66, 91, 113, 264, 175, 335, 109), 
            nrow = 32, ncol = 3, byrow = FALSE)

# Add column names
colnames(matrix_A) <- c("column_1", "column_2", "column_3")

使用 qr() 函数进行 QR 分解

我们可以用 qr() 函数一步完成该分解。将结果赋给 qr_decomp 对象,它是一个包含 QR 分解结果的列表。

# Perform and Display the QR decomposition
(qr_decomp <- qr(matrix_A))

当我们访问 QR 分解对象的秩(rank)组件时,可以看到分解矩阵中线性无关列的数量。输出表明该矩阵的秩为 3,意味着存在三列线性无关,这与我们的预期一致。访问 pivot 组件时,可以看到分解过程中使用的列置换信息。输出表明在分解过程中矩阵的列没有被置换,因为主元索引与原始列顺序一致。

qr_decomp$rank
qr_decomp$pivot
[1] 3
[1] 1 2 3

接下来,我们可分别用 qr.Q() 和 qr.R() 访问 Q 矩阵与 R 矩阵。输出显示 Q 矩阵为长方形,R 矩阵为方阵且对角线以下为零。

# Extract and Display the Q and R matrices
(Q_matrix <- qr.Q(qr_decomp))
(R_matrix <- qr.R(qr_decomp))
           [,1]         [,2]         [,3]
 [1,] -0.17721481  0.023462372 -0.015872336
 [2,] -0.17721481 -0.001889138 -0.060272525
 [3,] -0.19240465  0.077625112 -0.019039920
 [4,] -0.18059033 -0.030282834 -0.118470211
 [5,] -0.15780557 -0.089157964  0.064671266
 [6,] -0.15274229 -0.099258795 -0.187022354
 [7,] -0.12067485 -0.161573766  0.277715036
 [8,] -0.20590673  0.012764996 -0.275768785
 [9,] -0.19240465 -0.004891570 -0.156509401
[10,] -0.16202497 -0.082397566 -0.117343880
[11,] -0.15021065 -0.101326680 -0.120852866
[12,] -0.13839633 -0.182888939 -0.033166125
[13,] -0.14599125 -0.136918208  0.028289904
[14,] -0.12826977 -0.170282764  0.014320506
[15,] -0.08776353 -0.381326588 -0.165554550
[16,] -0.08776353 -0.398625266 -0.160607236
[17,] -0.12405037 -0.332631845 -0.083208412
[18,] -0.27341714  0.219354877 -0.069242665
[19,] -0.25653954  0.250472633 -0.021737724
[20,] -0.28607534  0.275923528 -0.005454217
[21,] -0.18143421  0.045632512 -0.033447865
[22,] -0.13080141 -0.140377926 -0.045388530
[23,] -0.12826977 -0.135983661 -0.031340393
[24,] -0.11223605 -0.201937338  0.228196652
[25,] -0.16202497 -0.122661730 -0.004593473
[26,] -0.23037925  0.176744505 -0.035884026
[27,] -0.21940881  0.138786843  0.013273146
[28,] -0.25653954  0.260613237  0.211010230
[29,] -0.13333305 -0.101525496  0.418085457
[30,] -0.16624437 -0.009027316  0.183837006
[31,] -0.12658201 -0.152109209  0.596664759
[32,] -0.18059033  0.012963861 -0.046253102
      column_1  column_2  column_3
[1,] -118.5003 -16.11602 -711.9200
[2,]    0.0000 -10.05857 -496.9319
[3,]    0.0000   0.00000  283.7369

最后,出于好奇,我们可以查看 qr() 函数用于分解的方法。

?qr

我们可以在输出中看到如下内容:

The function computes the QR decomposition of a matrix using the LINPACK or LAPACK routine dgeqrf.

如果我们查阅 FORTRAN 文档,会发现 dgeqrf 是一个使用 Householder 变换的 LAPACK 例程。

使用 Gram-Schmidt 过程进行 QR 分解

让我们尝试通过复现一个理念相近(非完全相同实现)的方法来理解 qr() 函数背后的原理,即按照 Gram-Schmidt 过程来思考。这个练习的目的不是精确复刻某种方法,而是帮助我们建立理解。

构造 Q 矩阵

为用 Gram-Schmidt 过程开始构造 Q 矩阵,我们需要先对子集数据的第一列进行归一化,得到 Q 的第一列。目标是将矩阵 A 中的一组向量转化为单位长度为一的新向量集合,并确保它们彼此正交。

要找到 Q 的第一列,我们提取 A 的第一列,计算欧几里得范数,然后进行归一化。

回顾一下,欧几里得范数衡量的是欧几里得空间中向量的长度或大小。我们把每个向量元素视为多维空间中的坐标,考虑从原点到向量端点的距离。对于向量 v = [v1, v2, …, vn],其欧几里得范数按以下方式计算:

欧几里得范数公式

因此,在 R 中实现时,我们对每一列求平方、求和、开方,最后进行归一化。

# Extract the first column of the subset data
col_1 <- matrix_A[, 1]

# Calculate the Euclidean norm of the first column
norm_col_1 <- sqrt(sum(col_1^2))

# Normalize and display the first column to get the first column of Q
(q1 <- col_1 / norm_col_1)

我们需要复习的第二个线性代数概念是投影。投影是一种将向量投影到子空间上的数学运算。子空间可被视为较大向量空间中的低维空间。投影的结果是一个位于子空间内、且距离原向量最近的向量,同时与该子空间正交。

我们用下图来帮助可视化,其中投影 p 是从 b 到 a 的最近点,并且必须与 b 垂直。

线性代数中的投影示意图

投影示意图。图片作者自制。

在 QR 分解的下一步中,我们首先计算矩阵 A 第二列在 Q 第一列方向上的投影,然后用 A 的第二列减去该投影,得到第二个与第一列正交的向量。最后对该向量归一化,得到 Q 的第二列。

我们来回顾一下投影的公式。投影矩阵等于 a 乘以 a 的转置,再除以标量 a 的转置乘以a,最后再乘以 b。

投影矩阵公式

或者:

投影矩阵公式

我们或许会问:这两种方法为何等价?当我们查看上面的第一个公式时,b 左侧的部分计算了一个矩阵,本质上捕捉了由向量 a 定义的方向及其缩放因子。这与在点积方法中计算 a 的单位向量是类似的。

下面的 R 代码使用点积方法计算向量 col_2 在向量 q1 方向上的投影。

# Extract the second column of the subset data
col_2 <- matrix_A[, 2]

# Calculate the projection of col_2 onto q1
projection_col_2_q1 <- sum(col_2 * q1) * q1

# Calculate the second orthogonal vector
v2 <- col_2 - projection_col_2_q1

# Calculate the Euclidean norm of the second orthogonal vector
norm_v2 <- sqrt(sum(v2^2))

# Normalize the second orthogonal vector to get the second column of Q
q2 <- v2 / norm_v2

如上代码所示,构造第二个正交向量以及 Q 矩阵对应的第二列,涉及投影、相减与归一化的过程。基本步骤与创建 Q 的第一列类似,但稍有不同且更为繁琐:我们先提取 A 的第二列;再对 Q 的第一列做投影;从 A 的第二列中减去该投影;接着进行归一化;最后得到 Q 的第二列。

我们可以把这一步对 q3 乃至最后一列重复。在我们的例子中,正交向量 q1、q2 和 q3 共同构成正交矩阵 Q 的基。这些向量很关键,因为它们为用一套新的坐标系来表达原始数据奠定了基础。

在这里我们看到,Gram-Schmidt 是一个迭代过程,确保每一个后续向量都与前面的向量正交。我们也要记住,正交意味着两个向量的点积为零,表示它们在向量空间中彼此垂直。我们用点积来验证 q1 与 q2 的正交性:

# Calculate the dot product of q1 and q2
dot_product_q1_q2 <- sum(q1 * q2)

# Print the dot product to confirm orthogonality
print(dot_product_q1_q2)
[1] 2.43295e-16

构造 R 矩阵

现在我们已经得到了 Q 矩阵,需要构造 R 矩阵。QR 分解中的 R 矩阵包含正交矩阵 Q 的各列与原始数据之间关系的信息。它是一个上三角矩阵,R 上三角部分的每个元素都由 Q 的相应列与原数据列之间的点积得到。

要构造上三角矩阵 R,我们需要计算 Q 的每一列与原始数据每一列的点积。由于之前的处理,最终得到的矩阵行数等于 Q 的列数,并呈上三角形式。

matrix_A <- matrix_A %>% as.matrix()

# Calculate and display the R matrix by multiplying Q^T and the subset data
R_matrix <- t(Q_matrix) %*% matrix_A

我们看到这次得到的 R 矩阵输出略有不同。原因在于 qr() 函数使用了数值稳定性更高的 Householder 反射,而非我们自己实现的、受 Gram-Schmidt 启发的过程。

         column_1      column_2  column_3
[1,] -1.185003e+02 -1.611602e+01 -711.9200
[2,]  5.814446e-16 -1.005857e+01 -496.9319
[3,]  7.633287e-15  1.338669e-15  283.7369

我们可以通过相乘 Q 与 R 并查看是否能还原 matrix_A 来验证重构。这里,%*% 是 R 中的矩阵乘法。

Q_matrix %*% R_matrix
   column_1 column_2 column_3
 [1,]     21.0    2.620      110
 [2,]     21.0    2.875      110
 [3,]     22.8    2.320       93
 [4,]     21.4    3.215      110
 [5,]     18.7    3.440      175
 [6,]     18.1    3.460      105
 [7,]     14.3    3.570      245
 [8,]     24.4    3.190       62
 [9,]     22.8    3.150       95
[10,]     19.2    3.440      123
[11,]     17.8    3.440      123
[12,]     16.4    4.070      180
[13,]     17.3    3.730      180
[14,]     15.2    3.780      180
[15,]     10.4    5.250      205
[16,]     10.4    5.424      215
[17,]     14.7    5.345      230
[18,]     32.4    2.200       66
[19,]     30.4    1.615       52
[20,]     33.9    1.835       65
[21,]     21.5    2.465       97
[22,]     15.5    3.520      150
[23,]     15.2    3.435      150
[24,]     13.3    3.840      245
[25,]     19.2    3.845      175
[26,]     27.3    1.935       66
[27,]     26.0    2.140       91
[28,]     30.4    1.513      113
[29,]     15.8    3.170      264
[30,]     19.7    2.770      175
[31,]     15.0    3.570      335
[32,]     21.4    2.780      109

线性回归中的 QR 分解

在 R 中,每次调用 lm() 函数都会使用 QR 分解。原因在于 QR 分解是一种数值稳定且高效的线性系统求解方法,尤其适用于大型数据集或当自变量(也就是设计矩阵的各列)近似共线时。

我们可以通过在控制台输入 lm 并查看函数签名中的 method = "qr" 和 qr = TRUE 来了解 lm() 如何使用 QR 分解。

lm
lm <- function (formula, data, subset, weights, na.action,
		method = "qr", model = TRUE, x = FALSE, y = FALSE,
		qr = TRUE, singular.ok = TRUE, contrasts = NULL,
		offset, ...)

还有其他方法可以求线性回归的系数。例如,我们可以通过求解所谓的正规方程来获得线性回归系数。我们甚至可以用微积分方法,对平方误差和对系数求偏导并令其为零。在简单线性回归(仅一个自变量与一个因变量)的情况下,我们还可以用相关系数 r 乘以 y 的标准差与 x 的标准差之比来得到回归线的斜率。

然而,在所有这些方法中,QR 分解可能是最优选择,因为它在稳定性与效率之间取得了平衡。尤其是正规方程 ATA x = ATB,计算代价非常高。因为为了解 x,我们写成 x = ((ATA)-1) ATB。其中需要对 ATA 求逆的部分,其时间复杂度为 O(n3),意味着当矩阵规模加倍时,求逆所需时间会增加到原来的八倍。这使其不适用于高维问题。

结语

QR 分解并非只存在于线性代数中。它在数据科学和机器学习中有着广泛用途。通过理解并应用 QR 分解,您可以提升处理复杂问题的能力,例如应对大型数据集与线性回归估计器。

此外,如果您对分解类方法感兴趣,值得一提的是,分解在数据科学中是一个更为通用的理念。比如,时间序列分解是将时间序列拆分为趋势-周期、季节性和余项,这些组件加总或相乘后可还原为原始序列。再如傅里叶变换,它将图像分解为频率成分,这在信号处理中大有用处。

当您继续探索数据科学中的高级主题时,会发现 QR 分解等概念与其他基础领域密切相关。欢迎学习我们的课程:R 中的数据科学线性代数 和 统计学入门。掌握这些相互关联的主题,将帮助您应对更具挑战性、也更丰富多样的数据科学问题。

Josef Waples

作者

Josef Waples

我是数据科学编辑,曾为科学期刊发表的研究论文撰稿。我的兴趣领域包括线性代数、统计学、R 等相关主题。

主题