请求改写RNimble的LQ分解函数以实现tensr包lq()功能

需求说明 #

• 目标：实现tensr包中lq()函数的功能，将m×n（m>n）矩形矩阵A分解为：
  • L：m×n下三角矩阵（主对角线以上元素为0）
  • Q：n×n正交矩阵（满足$Q^TQ=I_n$）
• 现有问题：当前代码对$A^T$做Gram-Schmidt分解时变量逻辑混乱，生成的是n×m正交矩阵和m×m下三角矩阵，不符合需求，需修正。

现有代码 #

nimQR.Q = nimbleFunction(
  run = function(A = double(2)){
    returnType(double(2))
    A1 = t(A)
    m <- dim(A1)[1]
    n <- dim(A1)[2]
    Q <- matrix(0, m, n)
    R <- matrix(0, n, n)
    for (k in 1:n) {
      Q[1:m, k] <- A1[1:m, k]
      if (k > 1) {
        for (i in 1:(k - 1)) {
          R[i, k] <- t(Q[1:m, i]) %*% Q[1:m, k]
          Q[1:m, k] <- Q[1:m, k] - R[i, k] * Q[1:m, i]
        }
      }
      R[k, k] <- sum(abs(Q[, k])^2)^(1/2)
      Q[1:m, k] <- Q[1:m, k]/R[k, k]
    }
    return(Q)
})

nimQR.R = nimbleFunction(
  run = function(A = double(2)){
    returnType(double(2))
    A1 = t(A)
    m <- dim(A1)[1]
    n <- dim(A1)[2]
    Q <- matrix(0, m, m)
    R <- matrix(0, n, n)
    for (k in 1:n) {
      Q[1:m, k] <- A1[1:m, k]
      if (k > 1) {
        for (i in 1:(k - 1)) {
          R[i, k] <- t(Q[1:m, i]) %*% Q[1:m, k]
          Q[1:m, k] <- Q[1:m, k] - R[i, k] * Q[1:m, i]
        }
      }
      R[k, k] <- sum(abs(Q[1:m, k])^2)^(1/2)
      Q[1:m, k] <- Q[1:m, k]/R[k, k]
    }
    return(t(R))
    
  })

改写后的解决方案 #

核心思路 #

LQ分解可通过对$A^T$做QR分解推导：

1. 对$A^{T$（n×m矩阵）执行QR分解，得到$A}T = Q_{QR} R_{QR}$，其中$Q_{QR}$是n×n正交矩阵，$R_{QR}$是n×m上三角矩阵。
2. 转置等式得到$A = R_{QR}^T Q_{QR}^T$，此时：
   • $L = R_{QR}^T$：m×n下三角矩阵（因$R_{QR}$是上三角，转置后为下三角）
   • $Q = Q_{QR}^T$：n×n正交矩阵（正交矩阵的转置仍为正交矩阵）

改写后的代码 #

# 返回LQ分解中的n×n正交矩阵Q
nimLQ.Q = nimbleFunction(
  run = function(A = double(2)){
    returnType(double(2))
    # 获取矩阵A的维度：m行n列（m>n）
    m <- dim(A)[1]
    n <- dim(A)[2]
    # 处理A的转置矩阵A^T（n×m）
    A_t <- t(A)
    # 初始化QR分解所需的矩阵
    Q_qr <- matrix(0, n, n)  # n×n正交矩阵
    R_qr <- matrix(0, n, m)  # n×m上三角矩阵
    
    # 对A^T的前n列执行Gram-Schmidt正交化（对应A的前n行）
    for (k in 1:n) {
      # 初始化当前列
      Q_qr[, k] <- A_t[, k]
      # 减去已正交化列的投影
      if (k > 1) {
        for (i in 1:(k-1)) {
          R_qr[i, k] <- t(Q_qr[, i]) %*% Q_qr[, k]
          Q_qr[, k] <- Q_qr[, k] - R_qr[i, k] * Q_qr[, i]
        }
      }
      # 归一化当前列
      norm_val <- sqrt(sum(Q_qr[, k]^2))
      R_qr[k, k] <- norm_val
      Q_qr[, k] <- Q_qr[, k] / norm_val
    }
    
    # LQ分解的Q是Q_qr的转置
    Q <- t(Q_qr)
    return(Q)
  })

# 返回LQ分解中的m×n下三角矩阵L
nimLQ.L = nimbleFunction(
  run = function(A = double(2)){
    returnType(double(2))
    m <- dim(A)[1]
    n <- dim(A)[2]
    A_t <- t(A)
    Q_qr <- matrix(0, n, n)
    R_qr <- matrix(0, n, m)
    
    for (k in 1:n) {
      Q_qr[, k] <- A_t[, k]
      if (k > 1) {
        for (i in 1:(k-1)) {
          R_qr[i, k] <- t(Q_qr[, i]) %*% Q_qr[, k]
          Q_qr[, k] <- Q_qr[, k] - R_qr[i, k] * Q_qr[, i]
        }
      }
      norm_val <- sqrt(sum(Q_qr[, k]^2))
      R_qr[k, k] <- norm_val
      Q_qr[, k] <- Q_qr[, k] / norm_val
    }
    
    # LQ分解的L是R_qr前n列的转置（对应A的n列）
    L <- t(R_qr[, 1:n])
    return(L)
  })

验证说明 #

调用函数后可通过以下方式验证正确性：

# 生成测试矩阵
set.seed(123)
A <- matrix(rnorm(10*5), nrow=10, ncol=5)  # m=10, n=5

# 执行LQ分解
L <- nimLQ.L(A)()
Q <- nimLQ.Q(A)()

# 验证A ≈ L %*% Q
all.equal(A, L %*% Q, tolerance=1e-6)

# 验证Q是正交矩阵：Q^T %*% Q ≈ I_n
all.equal(t(Q) %*% Q, diag(5), tolerance=1e-6)

# 验证L是下三角矩阵：主对角线以上元素为0
all.equal(L[upper.tri(L)], rep(0, sum(upper.tri(L))), tolerance=1e-6)

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