为用 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)技术支持
在 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 的第二列。
# 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)技术支持
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 启发的过程。
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 分解是一种基础的矩阵分解技术,广泛应用于数据科学和机器学习的诸多领域。它能将复杂矩阵拆解为更简单的部分,因此在求解线性代数问题与提升算法稳定性方面尤为重要。
无论您是在进行普通最小二乘线性回归,还是用 Python 或 PySpark构建推荐引擎,理解 QR 分解都能显著增强您处理大型数据集与复杂数学模型的能力。
什么是 QR 分解?
对矩阵进行因式分解(分解)是指将该矩阵表示为两个或多个矩阵的乘积。换言之,矩阵分解允许您把一个矩阵拆解为更简单的组件,这些组件在进行矩阵乘法后可还原为原矩阵。我们之所以重视矩阵分解,是因为它能简化许多线性代数问题,使之更易求解。
QR 分解是一种特定的矩阵分解方法,用于将矩阵分解为正交矩阵Q与上三角矩阵R的乘积。我们常见的等式为 A = 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,两部分都会用到它。使用 qr() 函数进行 QR 分解
我们可以用
qr()函数一步完成该分解。将结果赋给qr_decomp对象,它是一个包含 QR 分解结果的列表。当我们访问 QR 分解对象的秩(rank)组件时,可以看到分解矩阵中线性无关列的数量。输出表明该矩阵的秩为 3,意味着存在三列线性无关,这与我们的预期一致。访问 pivot 组件时,可以看到分解过程中使用的列置换信息。输出表明在分解过程中矩阵的列没有被置换,因为主元索引与原始列顺序一致。
接下来,我们可分别用
qr.Q()和qr.R()访问 Q 矩阵与 R 矩阵。输出显示 Q 矩阵为长方形,R 矩阵为方阵且对角线以下为零。最后,出于好奇,我们可以查看
qr()函数用于分解的方法。我们可以在输出中看到如下内容:
如果我们查阅 FORTRAN 文档,会发现
dgeqrf是一个使用 Householder 变换的 LAPACK 例程。使用 Gram-Schmidt 过程进行 QR 分解
让我们尝试通过复现一个理念相近(非完全相同实现)的方法来理解
qr()函数背后的原理,即按照 Gram-Schmidt 过程来思考。这个练习的目的不是精确复刻某种方法,而是帮助我们建立理解。构造 Q 矩阵
为用 Gram-Schmidt 过程开始构造 Q 矩阵,我们需要先对子集数据的第一列进行归一化,得到 Q 的第一列。目标是将矩阵 A 中的一组向量转化为单位长度为一的新向量集合,并确保它们彼此正交。
要找到 Q 的第一列,我们提取 A 的第一列,计算欧几里得范数,然后进行归一化。
回顾一下,欧几里得范数衡量的是欧几里得空间中向量的长度或大小。我们把每个向量元素视为多维空间中的坐标,考虑从原点到向量端点的距离。对于向量 v = [v1, v2, …, vn],其欧几里得范数按以下方式计算:
因此,在 R 中实现时,我们对每一列求平方、求和、开方,最后进行归一化。
我们需要复习的第二个线性代数概念是投影。投影是一种将向量投影到子空间上的数学运算。子空间可被视为较大向量空间中的低维空间。投影的结果是一个位于子空间内、且距离原向量最近的向量,同时与该子空间正交。
我们用下图来帮助可视化,其中投影 p 是从 b 到 a 的最近点,并且必须与 b 垂直。
投影示意图。图片作者自制。
在 QR 分解的下一步中,我们首先计算矩阵 A 第二列在 Q 第一列方向上的投影,然后用 A 的第二列减去该投影,得到第二个与第一列正交的向量。最后对该向量归一化,得到 Q 的第二列。
我们来回顾一下投影的公式。投影矩阵等于 a 乘以 a 的转置,再除以标量 a 的转置乘以a,最后再乘以 b。
或者:
我们或许会问:这两种方法为何等价?当我们查看上面的第一个公式时,b 左侧的部分计算了一个矩阵,本质上捕捉了由向量 a 定义的方向及其缩放因子。这与在点积方法中计算 a 的单位向量是类似的。
下面的 R 代码使用点积方法计算向量
col_2在向量q1方向上的投影。如上代码所示,构造第二个正交向量以及 Q 矩阵对应的第二列,涉及投影、相减与归一化的过程。基本步骤与创建 Q 的第一列类似,但稍有不同且更为繁琐:我们先提取 A 的第二列;再对 Q 的第一列做投影;从 A 的第二列中减去该投影;接着进行归一化;最后得到 Q 的第二列。
我们可以把这一步对
q3乃至最后一列重复。在我们的例子中,正交向量q1、q2和q3共同构成正交矩阵 Q 的基。这些向量很关键,因为它们为用一套新的坐标系来表达原始数据奠定了基础。在这里我们看到,Gram-Schmidt 是一个迭代过程,确保每一个后续向量都与前面的向量正交。我们也要记住,正交意味着两个向量的点积为零,表示它们在向量空间中彼此垂直。我们用点积来验证
q1与q2的正交性:构造 R 矩阵
现在我们已经得到了 Q 矩阵,需要构造 R 矩阵。QR 分解中的 R 矩阵包含正交矩阵 Q 的各列与原始数据之间关系的信息。它是一个上三角矩阵,R 上三角部分的每个元素都由 Q 的相应列与原数据列之间的点积得到。
要构造上三角矩阵 R,我们需要计算 Q 的每一列与原始数据每一列的点积。由于之前的处理,最终得到的矩阵行数等于 Q 的列数,并呈上三角形式。
我们看到这次得到的 R 矩阵输出略有不同。原因在于
qr()函数使用了数值稳定性更高的 Householder 反射,而非我们自己实现的、受 Gram-Schmidt 启发的过程。我们可以通过相乘 Q 与 R 并查看是否能还原
matrix_A来验证重构。这里,%*%是 R 中的矩阵乘法。