Hello, Currently, determinant(A) calculates the determinant of 'A' by factorizing A=LU and computing prod(diag(U)) [or the logarithm of the absolute value]. The factorization is done by LAPACK routine DGETRF, which gives a status code INFO, documented [1] as follows:
*> INFO is INTEGER *> = 0: successful exit *> < 0: if INFO = -i, the i-th argument had an illegal value *> > 0: if INFO = i, U(i,i) is exactly zero. The factorization *> has been completed, but the factor U is exactly *> singular, and division by zero will occur if it is used *> to solve a system of equations. Accordingly, when INF0>0, determinant(A) behaves as det(A)=0, _not_ computing prod(diag(U)). The problem here is that DGETRF can _also_ give positive INFO for matrices containing NaN, which may very well not be singular for some finite value of NaN. I claim that, when INFO>0, determinant(A) should _not_ behave as det(A)=0 unconditionally, but rather sometimes (depending on some test) give NaN. Here is one case where 0 is really "wrong": > (A <- diag(c(NaN, 1))) [,1] [,2] [1,] NaN 0 [2,] 0 1 > det(A) [1] 0 R isn't consistent, either: > (B <- diag(c(1, NaN))) [,1] [,2] [1,] 1 0 [2,] 0 NaN > det(B) [1] NaN Here, DGETRF _does_ succeed, because it does not "see" the trailing NaN in 'B'. So: Should R change to better handle the INFO>0 case? If so, how? Ideally (I think), the proposed change would give NaN for 'A' and 'B' above and 0 for 'C' and 'D' below (both of which really _are_ singular): > (C <- matrix(c(NaN, NaN, 0, 0), 2L, 2L)) [,1] [,2] [1,] NaN 0 [2,] NaN 0 > det(C) [1] NaN > (D <- t(C)) [,1] [,2] [1,] NaN NaN [2,] 0 0 > det(D) [1] 0 Furthermore, the proposed change should _not_ decrease the performance of determinant(A) for nonsingular 'A' ... For those looking, the relevant C-level function is det_ge_real(), defined in R-devel/src/modules/lapack/Lapack.c (at line 1260 in r83320). Mikael [1] https://github.com/Reference-LAPACK/lapack/blob/master/SRC/dgetrf.f ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel