gpt4 book ai didi

r - R 中的 chol() 函数不断返回上三角(我需要下三角)

转载 作者:行者123 更新时间:2023-12-01 10:32:53 25 4
gpt4 key购买 nike

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

最佳答案

我们可以使用 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?详情。

关于r - R 中的 chol() 函数不断返回上三角(我需要下三角),我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/40500998/

25 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com