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

[Mikael Jagan]

Hmm ... I can no longer reproduce this under r83996 when configuring
--without-lapack and --without-blas {the default}.  Now det(A), det(B),
det(C), and det(D) are all NaN.

I assume that the reason is the recent update to use LAPACK 3.11.0?
But I don't see any related NEWS here:

     https://netlib.org/lapack/lapack-3.11.0.html

It is possible that the underlying issue in DGETRF (INFO>0 due to NaN
pivots rather than 0 pivots {=> singular}) still exists for matrices
less trivial than my 2-by-2 examples.  Will have to experiment ...

Mikael

On 2022-11-09 3:58 pm, Mikael Jagan wrote:
> 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