[Rd] Proposal for diagonal-only chol2inv option

[Adam marstrand]

Dear R-devel,

I have a use-case where diag(chol2inv(U)) is required for many millions of
dense Cholesky factors. Currently, chol2inv() calls Lapack's DPOTRI, which
forms the complete inverse although only the diagonal may be needed.

After DTRTRI("U", "N", ...), each diagonal element of the inverse is the
squared 2-norm of one row of the inverted upper-triangular factor and can
therefore be obtained simply through using BLAS DDOT. This makes the work
after the triangular inversion quadratic instead of running the full DLAUUM
product.

It should be rather simple to implement, for example as chol2inv(x, size =
NCOL(x), diag.only = FALSE), with the relevant part of La_chol2inv()
becoming approximately:
if (only_diag) {
F77_CALL(dtrtri)("U", "N", &sz, REAL(ans), &sz, &info
       FCONE FCONE);
} else{
F77_CALL(dpotri("U", &sz, REAL(ans), &sz. &info FCONE);
}

if (only_diag) {
SEXP d = PROTECT(allocVector(REALSXP, sz)); nprot++;
BLAS_INT inc = (BLAS_INT) sz;
for (int i = 0; i < sz; i++) {
BLAS_INT len = (BLAS_INT) (sz - i);
double *x = REAL(ans) + i + SZ * i;
REAL(d)[i] = F77_CALL(ddot)(&len, x, &inc, x, &inc);
}
ans = d;
} else {
for (int j = 0; j < sz; j++)
for (int i = j + 1; i < sz; i++)
REAL(ans)[i + j * SZ] = REAL(ans)[j + i * SZ];
}

On my computer (current R-devel build/Windows/LAPACK 3.12.1), this cut the
duration of the operation almost in half (~1.8x faster) in my application
and in some quick benchmarking. It returns identical results to
diag(chol2inv(U)) from my testing.

Would such an option be of interest? I would be happy to provide a small
patch + test.

Best,
Adam Bøgh Marstrand-Jørgensen

	[[alternative HTML version deleted]]