R语言非满秩矩阵QR分解异常问题及正确实现需求

问题原因 #

原代码的偏差来自R默认qr()函数的列主元QR分解：分解结果满足的是 mat[, qr_decomp$pivot] = Q %*% R，而非直接的mat = Q %*% R。同时，若未利用秩信息和列置换，会导致对R矩阵维度的误解。

解决方法 #

方案一：列主元分解（数值稳定，推荐）

利用分解返回的列置换索引还原矩阵，确保误差小于1e-8：

set.seed(1)
n <- 30

X <- rnorm(n, mean = 0, sd = 1)
Z <- rnorm(n, mean = 0, sd = 1)
A <- rnorm(n, mean = 0, sd = 1)
W <- 4*X + 8*Z 
B <- 6*X + 9*Z + 7*A

mat <- cbind(X, Z, W, B, A)

# 执行列主元QR分解
qr_decomp <- qr(mat)
r <- qr_decomp$rank # 获取矩阵秩（结果为3）
Q <- qr.Q(qr_decomp) # 30×3的正交矩阵
R <- qr.R(qr_decomp) # 3×5的上三角矩阵（维度符合秩的要求）
pivot <- qr_decomp$pivot # 列置换的索引

# 按置换索引还原原矩阵
reconstructed_mat <- matrix(0, nrow = n, ncol = ncol(mat))
reconstructed_mat[, pivot] <- Q %*% R

# 验证误差
max_error <- max(abs(reconstructed_mat - mat))
cat("最大误差：", max_error, "\n") # 输出远小于1e-8

方案二：无列主元分解（仅适用于条件数较好的矩阵）

关闭列主元后，分解结果直接满足mat = Q %*% R：

set.seed(1)
n <- 30

X <- rnorm(n, mean = 0, sd = 1)
Z <- rnorm(n, mean = 0, sd = 1)
A <- rnorm(n, mean = 0, sd = 1)
W <- 4*X + 8*Z 
B <- 6*X + 9*Z + 7*A

mat <- cbind(X, Z, W, B, A)

# 关闭列主元的QR分解
qr_decomp <- qr(mat, pivot = FALSE)
r <- qr_decomp$rank
Q <- qr.Q(qr_decomp, complete = FALSE) # 30×3正交矩阵
R <- qr.R(qr_decomp) # 3×5上三角矩阵

# 验证误差
max_error <- max(abs(Q %*% R - mat))
cat("最大误差：", max_error, "\n") # 输出小于1e-8

说明 #

• 方案一的列主元分解数值稳定性更强，适合大多数实际场景；
• 两种方案得到的R矩阵维度均为3×5（与矩阵秩匹配），且还原后的矩阵与原矩阵误差远低于1e-8，满足数学上mat = Q %*% R的要求。

内容的提问来源于stack exchange，提问作者温泽海