如何在R中实现含固定列的矩阵分解？

已知矩阵y，需要将其分解为y = y1 %*% y2，其中y1的前k-1列固定为预设值，仅求解y1的最后一列和完整的y2。以下是具体的实现思路和代码：

问题数学定义 #

假设：

• y是n×m矩阵
• y1是n×k矩阵，前k-1列固定为矩阵A（n×(k-1)），最后一列为待求向量x（n×1），即y1 = cbind(A, x)
• y2是k×m矩阵，可拆分为(k-1)×m的矩阵B和m×1的向量c，即y2 = rbind(B, t(c))

分解式可展开为：
y = A %*% B + x %*% t(c)

我们的目标是求解x、B和c（即完整的y2）。

方法1：基于优化的最小二乘求解 #

当无法精确分解时，可通过最小化残差平方和获得近似解，使用R的optim()函数实现：

# 设定维度参数
n <- 5  # y的行数
m <- 3  # y的列数
k <- 3  # y1的列数（= y2的行数）

# 生成固定的y1前k-1列
A <- matrix(rnorm(n*(k-1)), nrow = n)

# 构造真实数据用于测试
x_true <- rnorm(n)
y2_true <- matrix(rnorm(k*m), nrow = k)
y <- A %*% y2_true[1:(k-1), ] + x_true %*% t(y2_true[k, ])

# 定义目标函数：计算残差平方和
objective <- function(params) {
  x <- params[1:n]
  y2 <- matrix(params[(n+1):(n + k*m)], nrow = k)
  y1 <- cbind(A, x)
  sum((y - y1 %*% y2)^2)
}

# 初始化参数
init_params <- c(rnorm(n), rnorm(k*m))

# 优化求解
result <- optim(init_params, objective, method = "L-BFGS-B")

# 提取结果
x_est <- result$par[1:n]
y2_est <- matrix(result$par[(n+1):(n + k*m)], nrow = k)

# 验证结果
y_est <- cbind(A, x_est) %*% y2_est
cat("优化法残差平方和：", sum((y - y_est)^2), "\n")

方法2：伪逆矩阵法（适用于近似/精确分解） #

利用矩阵伪逆先消除固定列A的影响，得到残差矩阵后提取共线分量：

library(MASS)

# 计算A的伪逆
A_pinv <- ginv(A)
# 计算残差矩阵（消除A的列空间影响）
R <- y - A %*% A_pinv %*% y

# 提取x的估计（残差矩阵的任意一列）
x_est2 <- R[, 1]
# 求解c：R = x_est2 %*% t(c)
c_est <- t(R) %*% x_est2 / sum(x_est2^2)
# 求解B：A%*%B = y - x_est2%*%t(c_est)
B_est <- A_pinv %*% (y - x_est2 %*% t(c_est))
# 组合得到y2
y2_est2 <- rbind(B_est, t(c_est))

# 验证结果
y_est2 <- cbind(A, x_est2) %*% y2_est2
cat("伪逆法残差平方和：", sum((y - y_est2)^2), "\n")

注意事项 #

• 如果是精确分解，两种方法都能得到精确解；如果是近似分解，优化法的结果更通用，伪逆法在残差矩阵秩为1时效率更高。
• 当矩阵规模较大时，可考虑使用更高效的优化算法或矩阵代数简化计算。

内容的提问来源于stack exchange，提问作者Nmgh