> The function ''det'' works improperly for a singular matrix
and returns a
> non-zero value even if ''solve'' reports singularity. The
matrix is very simple
> as shown below.
>
> A <- diag(rep(c(64,8), c(8,8)))
> A[9:16,1] <- 8
> A[1,9:16] <- 8
>
> det(A)
> #[1] -196608
> solve(A)
> #Error in solve.default(A) : system is computationally singular: reciprocal
> condition number = 2.31296e-18
The "det" function cannot work properly in the limited precision
of floating point numbers. May be, "det" could also do a test
of computational singularity and refuse to compute the value
similarly as "solve" does.
The eigen-values of your matrix are
> eigen(A)$values
[1] 7.200000e+01 6.400000e+01 6.400000e+01 6.400000e+01 6.400000e+01
[6] 6.400000e+01 6.400000e+01 6.400000e+01 8.000000e+00 8.000000e+00
[11] 8.000000e+00 8.000000e+00 8.000000e+00 8.000000e+00 8.000000e+00
[16] -2.023228e-15
The last value is not zero due to rounding. The determinant is the product
of eigenvalues, so we get something large.
The problem may also be seen as follows:
> det(A/8)
[1] -6.98492e-10
This is close to zero as it should be, however, det(A) = det(A/8)*8^16,
and this is what we really get:
> det(A/8)*8^16
[1] -196608
> det(A)
[1] -196608
There are even matrices closer to a diagonal matrix than A with
a similar problem:
> B <- 20*diag(16); B[1:2,1:2] <- c(25,35,35,49); det(B);
[1] 44565.24
All the best, Petr.