Daniel-Diaz / matrix

A Haskell native implementation of matrices and their operations.
BSD 3-Clause "New" or "Revised" License
35 stars 31 forks source link

inverse could check to see if matrix invertable #38

Open istathar opened 8 years ago

istathar commented 8 years ago

Given that the type signature of inverse is

inverse :: (Fractional a, Eq a) => Matrix a -> Either String (Matrix a)

I was a bit surprised that a matrix that isn't invertable didn't result in a Left:

λ> a
┌         ┐
│ 1.0 1.0 │
│ 2.0 2.0 │
└         ┘
λ> inverse a
Right ┌ *** Exception: getElem: Trying to get the (2,1) element from a 1x3 matrix.
CallStack (from HasCallStack):
  error, called at ./Data/Matrix.hs:434:16 in matrix-0.3.5.1-4UERfYC1lvKBLaK5oxpBHX:Data.Matrix
λ>

Boom.

λ> detLU a
0.0
λ> 

Oh.

I realize that large matrices don't suffer from this problem statistically, but smaller ones can bork pretty easily. So can we check? Calculating the determinant isn't free, but in the way you have a few different cases for matrix multiplication, can we have a safer inverse function that checks? You've already got the type signature to support it.

(I note in passing that there seem to be inverse methods that use the determinant in the calculation; Cayley–Hamilton and Cramer's Rule so if either of those were suitable for the implementation of inverse then the safety check calculating the determinant would contribute to the overall computation)

Anyway, I'd like inverse to not blow up if handed a non-invertable matrix. The other approach is to catch errors being thrown, but that's a bit ugly.

AfC

Daniel-Diaz commented 8 years ago

Thanks for the report.

Yes, I do think the inverse needs a review, or even a rewrite. It is not intentional that inverse is partial. It should always return a result. So this is a bug. However, I am not sure that calculating the determinant before calculating the inverse is a good approach. While calculating the inverse, it should become apparent at some step that the inverse can't be calculated, and therefore the input matrix was not invertible. If this is the case, the computation should stop there and return a Left value.

Daniel-Diaz commented 8 years ago

I added a test for this case that is currently failing.