R 中的 chol() 函数不断返回上三角（我需要下三角）答案

【问题标题】：chol() function in R keeps returning Upper Triangular (I need Lower Triangular)R 中的 chol() 函数不断返回上三角（我需要下三角）
【发布时间】：2016-11-09 06:14:26
【问题描述】：

我正在尝试使用chol() 函数在 R 中获取以下矩阵的下三角 Cholesky 分解。但是，它不断返回上三角分解，即使在查看文档之后，我似乎也找不到获得下三角分解的方法。以下是我正在使用的代码 -

x <- matrix(c(4,2,-2, 2,10,2, -2,2,5), ncol = 3, nrow = 3)
Q <- chol(x)
Q
#       [,1] [,2]      [,3]
#  [1,]    2    1 -1.000000
#  [2,]    0    3  1.000000
#  [3,]    0    0  1.732051

我基本上需要找到矩阵Q 这样QQ' = X。感谢您的帮助！

【问题讨论】：

但这会给我下三角分解吗？我试过这个，当我做 Q*t(Q) 它没有返回原始矩阵？
如果你先做L <- t(Q) 然后L %*% t(L) 你会得到原来的矩阵，这就是你想要的。
非常感谢！我刚刚意识到我使用 * 而不是 %*% 在 R 中乘以矩阵。哇。这是一个陡峭的学习曲线

标签： r matrix matrix-decomposition

【解决方案1】：

我们可以使用t()将得到的上三角转置得到下三角：

x <- matrix(c(4,2,-2, 2,10,2, -2,2,5), ncol = 3, nrow = 3)
R <- chol(x)  ## upper tri
L <- t(R)  ## lower tri
all.equal(crossprod(R), x)  ## t(R) %*% R
# [1] TRUE
all.equal(tcrossprod(L), x)  ## L %*% t(L)
# [1] TRUE

但是，我认为你不是唯一一个有这种疑问的人。因此，我将对此进行详细说明。

chol.default 来自 R 基础调用 LAPACK 例程 dpotrf 用于非旋转 Cholesky 分解，LAPACK 例程 dpstrf 用于旋转 Cholesky 分解。两个 LAPACK 例程都允许用户选择是使用上三角还是下三角，但是 R 禁用下三角选项并且只返回上三角。见源代码：

chol.default
#function (x, pivot = FALSE, LINPACK = FALSE, tol = -1, ...) 
#{
#    if (is.complex(x)) 
#        stop("complex matrices not permitted at present")
#    .Internal(La_chol(as.matrix(x), pivot, tol))
#}
#<bytecode: 0xb5694b8>
#<environment: namespace:base>

// from: R_<version>/src/modules/lapack.c
static SEXP La_chol(SEXP A, SEXP pivot, SEXP stol)
{
  // ...omitted part...
if(!piv) {
int info;
F77_CALL(dpotrf)("Upper", &m, REAL(ans), &m, &info);
if (info != 0) {
    if (info > 0)
    error(_("the leading minor of order %d is not positive definite"),
          info);
    error(_("argument %d of Lapack routine %s had invalid value"),
      -info, "dpotrf");
}
} else {
double tol = asReal(stol);
SEXP piv = PROTECT(allocVector(INTSXP, m));
int *ip = INTEGER(piv);
double *work = (double *) R_alloc(2 * (size_t)m, sizeof(double));
int rank, info;
F77_CALL(dpstrf)("U", &m, REAL(ans), &m, ip, &rank, &tol, work, &info);
if (info != 0) {
    if (info > 0)
    warning(_("the matrix is either rank-deficient or indefinite"));
    else
    error(_("argument %d of Lapack routine %s had invalid value"),
          -info, "dpstrf");
}
// ...omitted part...
return ans;
}

如您所见，它将“Upper”或“U”传递给 LAPACK。

上三角因子在统计学中更常见的原因是我们在最小二乘计算中经常将 Cholesky 因子分解与 QR 因子分解进行比较，而后者只返回上三角函数。要求 chol.default 始终返回上三角函数提供代码重用。例如，chol2inv 函数用于查找最小二乘估计的未缩放协方差，我们可以将chol.default 或qr.default 的结果提供给它。详情请见How to calculate variance of least squares estimator using QR decomposition in R?。

【讨论】：