2 of 3 edited body; edited tags; edited title

Type safe, easy to understand matrix multiplication

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.


module Data.Nat where

data Nat = Z | S Nat


  , 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)


  , 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
    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