Can good type systems distinguish between matrices in different bases?

Question

My program (Hartree-Fock/iterative SCF) has two matrices F and F' which are really the same matrix expressed in two different bases. I just lost three hours of debugging time because I accidentally used F' instead of F. In C++, the type-checker doesn't catch this kind of error because both variables are Eigen::Matrix<double, 2, 2> objects.

I was wondering, for the Haskell/ML/etc. people, whether if you were writing this program you would have constructed a type system where F and F' had different types? What would that look like? I'm basically trying to get an idea how I can outsource some logic errors onto the type checker.

Edit: The basis of a matrix is like the unit. You can say 1L or however many gallons, they both mean the same thing. Or, to give a vector example, you can say (0,1) in Cartesian coordinates or (1,pi/2) in polar. But even though the meaning is the same, the numerical values are different.

Edit: Maybe units was the wrong analogy. I'm not looking for some kind of record type where I can specify that the first field will be litres and the second gallons, but rather a way to say that this matrix as a whole, is defined in terms of some other matrix (the basis), where the basis could be any matrix of the same dimensions. E.g., the constructor would look something like mkMatrix [[1, 2], [3, 4]] [[5, 6], [7, 8]] and then adding that object to another matrix would type-check only if both objects had the same matrix as their second parameters. Does that make sense?

Edit: definition on Wikipedia, worked examples

Answer 1

This is entirely possible in Haskell.

Statically checked dimensions

Haskell has arrays with statically checked dimensions, where the dimensions can be manipulated and checked statically, preventing indexing into the wrong dimension. Some examples:

This will only work on 2-D arrays:

multiplyMM :: Array DIM2 Double -> Array DIM2 Double -> Array DIM2 Double 

An example from repa should give you a sense. Here, taking a diagonal requires a 2D array, returns a 1D array of the same type.

diagonal :: Array DIM2 e -> Array DIM1 e 

or, from Matt sottile's repa tutorial, statically checked dimensions on a 3D matrix transform:

f :: Array DIM3 Double -> Array DIM2 Double f u =   let slabX = (Z:.All:.All:.(0::Int))       slabY = (Z:.All:.All:.(1::Int))       u' = (slice u slabX) * (slice u slabX) +            (slice u slabY) * (slice u slabY)   in     R.map sqrt u' 

Statically checked units

Another example from outside of matrix programming: statically checked units of dimension, making it a type error to confuse e.g. feet and meters, without doing the conversion.

 Prelude> 3 *~ foot + 1 *~ metre  1.9144 m 

or for a whole suite of SI units and quanities.

E.g. can't add things of different dimension, such as volumes and lengths:

> 1 *~ centi litre + 2 *~ inch  Error: Expected type: Unit DVolume a1   Actual type: Unit DLength a0 

So, following the repa-style array dimension types, I'd suggest adding a Base phantom type parameter to your array type, and using that to distinguish between bases. In Haskell, the index Dim type argument gives the rank of the array (i.e. its shape), and you could do similarly.

Or, if by base you mean some dimension on the units, using dimensional types.

So, yep, this is almost a commodity technique in Haskell now, and there's some examples of designing with types like this to help you get started.

Answer 2

This is a very good question. I don't think you can encode the notion of a basis in most type systems, because essentially anything that the type checker does needs to be able to terminate, and making judgments about whether two real-valued vectors are equal is too difficult. You could have (2 v_1) + (2 v_2) or 2 (v_1 + v_2), for example. There are some languages which use dependent types [ wikipedia ], but these are relatively academic.

I think most of your debugging pain would be alleviated if you simply encoded the bases in which you matrix works along with the matrix. For example,

newtype Matrix = Matrix { transform :: [[Double]],      srcbasis :: [Double], dstbasis :: [Double] } 

and then, when you M from basis a to b with N, check that N is from b to c, and return a matrix with basis a to c.

NOTE -- it seems most people here have programming instead of math background, so I'll provide short explanation here. Matrices are encodings of linear transformations between vector spaces. For example, if you're encoding a rotation by 45 degrees in R^2 (2-dimensional reals), then the standard way of encoding this in a matrix is saying that the standard basis vector e_1, written "[1, 0]", is sent to a combination of e_1 and e_2, namely [1/sqrt(2), 1/sqrt(2)]. The point is that you can encode the same rotation by saying where different vectors go, for example, you could say where you're sending [1,1] and [1,-1] instead of e_1=[1,0] and e_2=[0,1], and this would have a different matrix representation.

Edit 1

If you have a finite set of bases you are working with, you can do it...

{-# LANGUAGE EmptyDataDecls #-} data BasisA data BasisB data BasisC  newtype Matrix a b = Matrix { coefficients :: [[Double]] }  multiply :: Matrix a b -> Matrix b c -> Matrix a c multiply (Matrix a_coeff) (Matrix b_coeff) = (Matrix multiplied) :: Matrix a c     where multiplied = undefined -- your algorithm here 

Then, in ghci (the interactive Haskell interpreter),

*Matrix> let m = Matrix [[1, 2], [3, 4]] :: Matrix BasisA BasisB *Matrix> m `multiply` m  <interactive>:1:13:     Couldn't match expected type `BasisB'         against inferred type `BasisA' *Matrix> let m2 = Matrix [[1, 2], [3, 4]] :: Matrix BasisB BasisC *Matrix> m `multiply` m2 -- works after you finish defining show and the multiplication algorithm