[Rd] det(diag(c(NaN, 1))) should be NaN, not 0

Mikael Jagan j@g@nmn2 @end|ng |rom gm@||@com
Wed Nov 9 21:58:52 CET 2022


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



More information about the R-devel mailing list