> 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. ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
