Haskell Matrices testing using QuickCheck

Question

I'm creating a Matrix module in Haskell and I want to use QuickCheck to test some properties of my code. Specifically I want to generate random matrices that have an associated inverse. The following is my attempt at creating a QuickCheck generator that generates such matrices.

invertibleMatrix :: (Num a, Arbitrary a) => Gen (Matrix a)
invertibleMatrix = do s <- choose (2,10)
                      a <- vectorOf s (vector s)
                      if (det (Matrix a) == 0) then
                        invertibleMatrix
                      else
                        return (Matrix a)

The code first creates a size between 2 and 10 and then a vector of vectors of this size. If the determinant is zero then the matrix is not invertible and so I call invertibleMatrix recursively. Otherwise I return the new matrix.

The trouble is, this code does not terminate when I put it in a property to test. (I think that it's constantly creating the same s x s matrix of zero elements which obviously doesn't have an inverse and so it goes into an infinite loop). What am I doing wrong? How do I fix this? Thanks.

Mark

Answer 1

squareMatrix :: (Num a, Arbitrary a) => Gen (Matrix a)
squareMatrix = do s <- choose (2,6)
                  a <- vectorOf s (vector s)
                  return (Matrix a)

invertibleMatrix :: (Num a, Arbitrary a) => Gen (Matrix a)
invertibleMatrix = suchThat squareMatrix ((/=0) . det) 

In case anyone wants to know, this is the solution.

Mark

Answer 2

As a way of bypassing your problem, you could notice that if A is an n×n matrix then A - tI is invertible for all but at most n values of t (i.e. the eigenvalues of A). So generate a matrix, if it is not invertible add a small multiple of the identity to it, and keep trying until it is. Then the process is guaranteed to terminate (as long as the underlying numeric type behaves reasonably well, which floating points sometimes won't, e.g. if the entries of A are much larger than the values of t that you try).

Answer 3

I have not tested this, but I think that the sort of perturbation you want can be taken from CoArbitrary which is usually used to implement a random function.

invertibleMatrix :: (Num a, Arbitrary a) => Gen (Matrix a)
invertibleMatrix = kick seed where
  seed :: Integer
  seed = 0
  kick n = do
    s <- choose (2,10)
    a <- vectorOf s (vector s)
    if (det (Matrix a) == 0) then
      coarbitrary (kick (succ n))
    else
      return (Matrix a)