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

Question

Я пытаюсь получить разложение по нижнему треугольнику Холецского следующей матрицы в 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 для неповоротной факторизации Холецкого и рутины 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.

Причина, по которой верхний треугольный фактор чаще встречается в статистике, состоит в том, что мы часто сравниваем факторизацию Холецкого с QR-факторизацией в вычислении методом наименьших квадратов, в то время как последний возвращает только верхний треугольник. требующий chol.default всегда возвращать верхний треугольник предлагает повторное использование кода. Например, chol2inv Функция используется для нахождения немасштабированной ковариации оценки наименьших квадратов, где мы можем скормить ее либо результатом chol.default или же qr.default, См. Как рассчитать дисперсию оценки наименьших квадратов, используя QR-разложение в R? для деталей.