Функция chol() в R продолжает возвращать верхний треугольник (мне нужен нижний треугольник)

Question

Я пытаюсь получить нижестоящий треугольный 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. Спасибо за вашу помощь!

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 базовых вызовов 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; 
} 

Как вы можете видеть, она проходит «верхний» или «U» в LAPACK.

Причина, по которой верхний треугольный фактор чаще встречается в статистике, заключается в том, что мы часто сравниваем факторизацию Холецкого с QR-факторизацией при вычислении наименьших квадратов, в то время как последний возвращает верхний треугольник. Требование chol.default всегда возвращать верхние треугольные предложения повторного использования кода. Например, функция chol2inv используется для поиска немасштабированной ковариации наименьшей квадратной оценки, где мы можем ее подавать либо в результате chol.default, либо qr.default. См. How to calculate variance of least squares estimator using QR decomposition in R?.