R中的chol（）函数不断返回上三角形（我需要下三角形）

我试图使用chol()函数在R中得到下面的矩阵的下三角Cholesky分解。然而，它不断返回上三角分解，我似乎无法找到一种方法来获得下三角分解，即使在查看文档之后。下面是我使用的代码 -R中的chol（）函数不断返回上三角形（我需要下三角形）

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。谢谢你的帮助！

来源

2016-11-09 tattybojangler

但会给我下三角分解？我试过这个，当我做了Q * t（Q）它没有返回原始矩阵？ – tattybojangler

如果你做'L < - t（Q）'，然后'L％*％t（L）'，你可以得到你想要的原始矩阵。 – Bhas

非常感谢！我刚刚意识到我使用*而不是％*％来在R.Wow中乘以矩阵。这是一个陡峭的学习曲线 – tattybojangler

我们可以使用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的结果。详情请参阅How to calculate variance of least squares estimator using QR decomposition in R?。

来源

2016-11-09 06:28:08

R中的chol（）函数不断返回上三角形（我需要下三角形）

回答

相关问题