Type-safe matrix multiplication

StackOverflow https://stackoverflow.com/questions/8332392

Question

After the long-winded discussion at Write this Scala Matrix multiplication in Haskell, I was left wondering...what would a type-safe matrix multiplication look like? So here's your challenge: either link to a Haskell implementation, or implement yourself, the following:

data Matrix ... = ...

matrixMult :: Matrix ... -> Matrix ... -> Matrix ...
matrixMult ... = ...

Where matrixMult produces a type error at compile time if you try to multiply two matricies with incompatible dimensions. Brownie points if you link to papers or books that discuss this precise topic, and/or discuss yourself how useful/useless this functionality is.

Was it helpful?

Solution

There are a number of packages that implement this:

The Repa papers in particular have a really nice discussion of the design space and choices made: http://repa.ouroborus.net/

Of historical interest is McBride's "Faking It" from 2001 which describes strongly typed vectors. The techniques he employs are fairly similar to those used in the above packages. They were obviously known in circles doing dependently typed programming, but my impression is that the "Faking It" paper is one of the earlier instances where these were used in Haskell. Oleg's 2005 Monad Reader article on number-parameterized types has some good discussion on the history of these techniques as well.

OTHER TIPS

You could use type-level natural numbers to encode the dimensions. Your matrix type becomes

-- x, y: Dimensions
data Matrix x y n = ...

and you have to define two additonal ADTs and a class TLN (Type Level Naturals):

data Zero
data Succ a
class    TLN a                 where fromTLN :: a -> Int
instance TLN Zero              where fromTLN = const Zero
instance TLN a => TLN (Succ a) where fromTLN = 1 + fromTLN (undefined :: a)

Your function's type is quite easy:

matrixMult :: (TLN x, TLN y, TLN t, Num a) =>
  Matrix x t a -> Matrix t y a -> Matrix x y a

You can extract the arrays dimension by generating an undefined of appropriate type together with the ScopedTypeVariables extension.

This code is completely untested and GHC may barf on compilation. It is just a sketch about how it could be done.

Link

Sorry, can't resist pasting something I whipped up ages ago. This was before type families so I used fundeps for the arithmetic. I verified that this still works on GHC 7.

{-# LANGUAGE EmptyDataDecls,
  ScopedTypeVariables,
  MultiParamTypeClasses,
  FunctionalDependencies,
  FlexibleContexts,
  FlexibleInstances,
  UndecidableInstances #-}

import System.IO


-- Peano type numerals

data Z
data S a

type One = S Z
type Two = S One
type Three = S Two

class Nat a
instance Nat Z
instance Nat a => Nat (S a)

class Positive a
instance Nat a => Positive (S a)

class Pred a b | a -> b
instance Pred (S a) a


-- Vector type

newtype Vector n k = Vector {unVector :: [k]}
    deriving (Read, Show, Eq)

empty :: Vector Z k
empty = Vector []

vtail :: Pred s' s => Vector s' k -> Vector s k
vtail (Vector (a:as)) = Vector as
vhead :: Positive s => Vector s k -> k
vhead (Vector (a:as)) = a

liftV :: (a->b) -> Vector s a -> Vector s b
liftV f = Vector . map f . unVector

type Matrix m n k = Vector m (Vector n k)

infixr 6 |>
(|>) :: k -> Vector s k -> Vector (S s) k
k |> v = Vector . (k:) . unVector $ v


-- Arithmetic

instance (Num k) => Num (Vector n k) where
    (+) (Vector v) (Vector u) = Vector $ zipWith (+) v u
    (*) (Vector v) (Vector u) = Vector $ zipWith (*) v u
    abs = liftV abs
    signum = liftV signum

dot :: Num k => Vector n k -> Vector n k -> k
dot u v = sum . unVector $ v*u

class Transpose n m where
    transpose :: Matrix n m k -> Matrix m n k

instance (Transpose m a, Nat a, Nat m) => Transpose m (S a) where
    transpose v = liftV vhead v |>
                  transpose (liftV vtail v)

instance Transpose m Z where
    transpose v = empty

multiply :: (Nat n, Nat m, Nat n', Num k, Transpose m n) =>
            Matrix m n k -> Matrix n' m k -> Matrix n n' k
multiply a (Vector bs) = Vector [Vector [a `dot` b | a <- as] | b <- bs]
    where (Vector as) = transpose a

printMatrix :: Show k => Matrix m n k -> IO ()
printMatrix = mapM_ (putStrLn) . map (show.unVector) . unVector


-- Examples

m :: Matrix Three Three Integer
m =    (1 |> 2 |> 3 |> empty)
    |> (2 |> 3 |> 4 |> empty) 
    |> (3 |> 4 |> 5 |> empty) |> empty
n :: Matrix Three Two Integer
n =    (1 |> 0 |> empty)
    |> (0 |> 1 |> empty) 
    |> (1 |> 1 |> empty) |> empty
o = multiply n m
p = multiply n (transpose n)

It is more Haskell-idiomatic to speak not actually of matrices, whose dimension is just a number that tells you little about the structure of the mapping / the spaces it maps between. Instead, matrix multiplication is best treated as the category-composition operation in the category Vectk. Vector spaces are naturally represented by Haskell types; the vector-space library has this for a long time.

As a composition of linear functions, the dimensions-check is then a corollary of the type-checking that's done anyway for Haskell function compositions. And not only that, you can also distinguish between different spaces that may be incompatible despite having the same dimension – for instance, matrices themselves form a vector space (a tensor space), but the space of 3×3 matrices isn't really compatible with the space of 9-element vectors. In Matlab and other “array languages”, dealing with linear mappings on a space of linear mapping requires error-prone reshaping between tensors of different rank; surely we don't want this in Haskell!

There's one catch: to implement these functions efficiently, you can't just have functions between any spaces, but need a kind of underlying representation that's still matrix-like. That only works when all permitted spaces are actually vector spaces, so you can't use the standard Category class as that demands an id between any two types. Instead you need a category class that's actually restricted to vector spaces. That's not really hard to express in modern Haskell though.

Two libraries that have gone this route are:

  • Mike Izbicki's subhask, which uses hmatrix matrices internally but exposes a nice high-level interface.
  • My own linearmap-category, which uses a dedicated implementation of the tensors spanned by each space.
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top