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.

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
• Why no one wants to review this. =( – Ignat Insarov Mar 9 '18 at 3:45
• More code = more to read = more to review = needs more time. Also, you're using features that are seldom reviewed. There are 7 questions that use DataKinds, and only 3 of them have been reviewed yet, so it might take some time. That being said, what was your rational to use a newtype instead of a GADT for Vector? – Zeta Mar 9 '18 at 6:45
• @Zeta I do not understand the hype about generalized algebraic types, I do not see the benefits they supposedly bring, I did not want to blindly copycat advanced code. In other words, there is simply no rationale. // Looks like there's only a handful of people doing reviews of Haskell code here. – Ignat Insarov Mar 9 '18 at 9:50

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.

• I understand it that, these points notwithstanding, you find the code useable? Did I succeed at making matrices rectangular (not ragged) by definition and making matrix multiplication work only on matrices of compatible size? Are these devices (Vector, Matrix, apply, compose, etc) consistent with each other? If you could add some positive points to your review, even very briefly, it will reassure me much. – Ignat Insarov Mar 10 '18 at 15:24
• @IgnatInsarov for ease of use, smart constructors are missing. Otherwise the functions are consistent and the library seems sane. Note that compose is usually not implemented like this, unless you really want to multiply two matrices (which you try to avert as much as possible due to the $\mathcal O(n^3)$ complexity). Please note that this a "dry" review. I don't have GHC at hand at the moment, only repl.it. – Zeta Mar 11 '18 at 8:56
• There is a typo in a formula in the compose section. $A$ is defined two times. – Ignat Insarov Mar 11 '18 at 10:34
• Thanks, fixed. One last remark on ease of use: you mix examples and library functions in your code. If you just keep the library functions and add the example as additional codeblock/module, a follow-up question might get more traction. – Zeta Mar 11 '18 at 10:51