R中的chol()函数继续返回上三角(我需要下三角)
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。谢谢你的帮忙!
回答 1
Stack Overflow用户
回答已采纳
发布于 2016-11-09 06:28:08
我们可以使用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但是，我认为你并不是唯一一个有这样疑问的人。因此，我将更详细地阐述这一点。
R基的chol.default调用用于非旋转Cholesky因式分解的LAPACK例程dpotrf和用于旋转Cholesky分解的LAPACK例程dpstrf。两个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?。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/40500998
复制相关文章
相似问题
腾讯云开发者
