R中的chol()函数保持返回Upper Triangular(我需要Lower Triangular)

Question

我试图使用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.谢谢你的帮助!

Answer 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 base调用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.请参阅如何使用R中的QR分解计算最小二乘估计的方差？详情.