Type safe, easy to understand matrix multiplication

Question

I always hated matrices (as they are taught in those linear algebra course books that I have browsed), because it's an obscure way of saying "a tuple of linear functions of same arity". It is very hard to remember which matrices you can multiply, what side the input is, what side the output, what transpositions should be done. Now, (square) matrices are a major example of a group, and matrices as an encoding of linear transformations are indispensable in differential geometry of any kind − for just two examples; so I had to stand up to it and design the matrices I would be able to like.

Now, one straightforward way of representing matrices in Haskell would be a nested list. Unfortunately, dealing in nested lists has proven to multiply my hate, rather than my matrices. I had to come up with a safer approach. What will follow is a take on it.

I post this code because it is actually above my head: I don't understand dependent typing at all, and half the language extensions I had to use I switched on because GHC told me to, not because I understand why they are indispensable here. I'd like to know if I'm on the right track − whether this code is sensible or not. I am happy about it now; I'd like to become less happy.

Is this implementation actually as type safe as I think?

Data.Nat

module Data.Nat where

data Nat = Z | S Nat

Data.Vector

{-# LANGUAGE
    DataKinds
  , KindSignatures
  , FlexibleInstances
  , ScopedTypeVariables
  , MultiParamTypeClasses
  , FunctionalDependencies
  , UndecidableInstances
  #-}

module Data.Vector where

import Data.Nat

newtype Vector (n :: Nat) a = Vector { unVector :: [a] } deriving Show

-- `vector` is thanks to HTNV @ stackoverflow.com/a/49120144/2108477
class MkVector n a n' r | r -> n a n' where
  mkVector :: (Vector n a -> Vector n' a) -> r

instance MkVector Z a n (Vector n a) where
  mkVector f = f $ Vector []

instance (MkVector n a n' r) => MkVector (S n) a n' (a -> r) where
  mkVector f x = mkVector $ \(Vector xs) -> f $ Vector $ x:xs

vector :: MkVector n a n r => r
vector = mkVector id

v0 :: Vector Z a
v0 = Vector [ ]

put :: a -> Vector n a -> Vector (S n) a
put x (Vector xs) = Vector (x:xs)

instance Functor (Vector n) where
    fmap f (Vector xs) = Vector (fmap f xs)

instance Foldable (Vector n) where
    foldMap f (Vector xs) = foldMap f xs

instance Repeat (Vector n) => Applicative (Vector n) where
    pure x = rep x
    (Vector fs) <*> (Vector xs) = Vector (zipWith ($) fs xs)

class Repeat a where
    rep :: forall x. x -> a x

instance Repeat (Vector Z) where
    rep _ = v0

instance Repeat (Vector n) => Repeat (Vector (S n)) where
    rep (x :: x) = Vector $ x: unVector (rep x :: Vector n x)

Data.Matrix

{-# LANGUAGE
    DataKinds
  , KindSignatures
  , FlexibleInstances
  , FlexibleContexts
  , UndecidableInstances
  , ExistentialQuantification
  , StandaloneDeriving
  , DeriveFunctor
  #-}

module Data.Matrix where

import Data.Nat
import Data.Vector
import Data.List (foldl1')

-- | A linear function of n variables.
newtype LinearFunction (n :: Nat) a = LinearFunction (Vector n a)
    deriving (Show, Functor)

instance Repeat (Vector n) => Applicative (LinearFunction n) where
    pure x = LinearFunction $ pure x
    (LinearFunction fs) <*> (LinearFunction xs) = LinearFunction $ fs <*> xs

-- | A few functions of the same arity.
data Matrix inp outp a = Matrix (Vector outp (LinearFunction inp a))
                       | forall inter. Sequence (Matrix inp inter a) (Matrix inter outp a)

deriving instance Show a => Show (Matrix inp outp a)

v1 = vector 1 2
v2 = vector 3 5
m1 :: Matrix (S (S Z)) (S (S Z)) Integer
m1 = Matrix (vector (LinearFunction v1) (LinearFunction v2))

-- |
--
-- λ apply v1 (Sequence m1 (Sequence m1 m1))
-- Vector [191,493]
--
apply :: Num a => Vector inp a -> Matrix inp outp a -> Vector outp a
apply v (Matrix (Vector fs)) = Vector $ applyFunction v <$> fs
apply v (Sequence m1 m2) = v `apply` m1 `apply` m2

-- |
--
-- λ applyFunction v1 (LinearFunction v2)
-- 13
--
applyFunction :: Num a => Vector n a -> LinearFunction n a -> a
applyFunction (Vector xs) (LinearFunction (Vector qs)) = sum $ zipWith (*) xs qs

compose_trivial, compose
    :: (Num a, Repeat (Vector inter), Repeat (Vector inp))
    => Matrix inp inter a -> Matrix inter outp a -> Matrix inp outp a

compose_trivial = Sequence

-- |
--
-- λ v1 `apply` (m1 `compose` m1 `compose` m1)
-- Vector {unVector = [191,493]}
--
compose (Matrix fs) (Matrix gs) = Matrix $ fmap (createRow fs) gs
  where
    createRow fs (LinearFunction qs)
        = let (Vector xs) = (\q -> fmap (* q)) <$> qs <*> fs
          in  foldl1' (\xs ys -> (+) <$> xs <*> ys) xs
compose m1 m2@Sequence{} = Sequence m1 m2
compose m1@Sequence{} m2 = Sequence m1 m2

Answer 1

Module names

Your module names clash with modules in popular packages. Data.Vector is provided by vector, Data.Matrix by matrix, Data.Nat by data-nat, although the last isn't often used. If you ever want to use Data.Vector in your project (for constant indexing, for example) you will have to pull some tricks.

Instead either use some other names, or use Ignat or Insarov in your module path as long as you're not sure what final name you want to use, e.g.

module Ignat.Data.Vector

or

module Data.Ignat.Vector

or

module Data.Vector.Ignat

Type signatures

There are several top-level type bindings that are missing a signature.

Type safety

Is this implementation actually as type safe as I think?

It's possible to create Vector's with wrong sizes:

main = print (Vector [] :: Vector (S Z) Int)

That's because the n parameter in your newtype is a phantom. Therefore, we can even use coerce from Data.Coerce to change vectors into arbitrary other ones:

main = print (coerce (v0 :: Vector Z Int) :: Vector (S (S Z)) Int)

Now you might say that this is an obvious malicous attempt to break type safety. However, we don't need coerce to do that. We can just write our own variant by accident, for example to use some list functions:

listFunc :: [Int] -> [Int] -- some library function you want to use

broken :: Vector n a -> Vector k a
broken = Vector . listFunc . unVector

We can use listFunc = id or listFunc = take 2 to see that it's possible to convert between vectors arbitrarily. This means that we can, by a small unfortunate accident, create a function that will happily compile:

apply :: Num a => Vector inp a -> Matrix inp outp a -> Vector out a

Note the missing p in Vector outp a.

In order to make those modules type safe, you have to constrain the exports so that no Vector or Matrix can be build or matched with their constructors our of their respective module. You had a similar loophole in your prime code, by the way.

A Vector based on GADTs does not have this problem:

data Vector (n :: Nat) a where
  VNil  :: Vector Z a
  VCons :: a -> Vector n a -> Vector (S n) a

Fixed-Length Vector Types in Haskell by Justin Le is worth a read at that point.

compose

compose has a worst- case \$\mathcal O(nkm)\$ complexity due to the matrix-matrix multiplication. If possible, you never actually want to use matrix-matrix multiplication, only matrix-vector ones, since the asymptotic time complexity is a lot better. If we have

$$ A B x $$ with \$A \in \mathbb K^{n\times k} \$, \$B \in \mathbb K^{k\times m} \$, \$x \in \mathbb K^{m} \$, we end up with

$$ \mathcal O(apply(mult(AB),x)) = \mathcal O(nkm+nm) = \mathcal O(nkm), $$ a cubic algorithm, dominated by the matrix-matrix multiplication.

On the other hand, if we use compose_trivial, we end up with $$ \mathcal O(seq(mult(AB),x)) = \mathcal O(A(Bx)) = \mathcal O(nk + km), $$ two matrix-vector multiplications. For small n, k or m that's not a huge difference, but if you start to use your type-safe matrices and vectors for numerical simulation, we end up in trouble.

Also, matrix-matrix multiplication can create a dense matrix from two sparse ones, but that's going to far at this point.

By the way, compose_trivial does not need Repeat.

Overall

I like the general approach. The ideas hold, but for ease of use we need a way to create matrices and vectors easily:

vector :: [a]   -> Maybe (Vector (n :: Nat) a)
matrix :: [[a]] -> Maybe (Matrix (outer :: Nat) (inner :: Nat) a)

Now add export lists and you have a nice, type safe vector/matrix multiplication library. However, I'm not sure whether S (S (S (… (S Z) …))) for large matrices leads to problems in GHC.