Write this Scala Matrix multiplication in Haskell

Question

Possible Duplicate:
Can you overload + in haskell?

Can you implement a Matrix class and an * operator that will work on two matrices?:

scala> val x = Matrix(3, 1,2,3,4,5,6)  
x: Matrix =   
[1.0, 2.0, 3.0]  
[4.0, 5.0, 6.0]  

scala> x*x.transpose  
res0: Matrix =   
[14.0, 32.0]  
[32.0, 77.0]  

and just so people don't say that it's hard, here is the Scala implementation (courtesy of Jonathan Merritt):

class Matrix(els: List[List[Double]]) {  

    /** elements of the matrix, stored as a list of  
      its rows */  
    val elements: List[List[Double]] = els  
    def nRows: Int = elements.length  
    def nCols: Int = if (elements.isEmpty) 0   
                     else elements.head.length  

    /** all rows of the matrix must have the same  
        number of columns */  
    require(elements.forall(_.length == nCols))  

    /* Add to each elem of matrix */
    private def addRows(a: List[Double],   
                        b: List[Double]):   
      List[Double] =   
        List.map2(a,b)(_+_)  

    private def subRows(a: List[Double],   
                        b: List[Double]):List[Double] =  
        List.map2(a,b)(_-_)  

    def +(other: Matrix): Matrix = {  
      require((other.nRows == nRows) &&   
              (other.nCols == nCols))  
      new Matrix(  
        List.map2(elements, other.elements)  
          (addRows(_,_))
      )  
    }  

    def -(other: Matrix): Matrix = {  
      require((other.nRows == nRows) &&   
              (other.nCols == nCols))  
      new Matrix(  
        List.map2(elements, other.elements)  
          (subRows(_,_))  
      )  
    }  

    def transpose(): Matrix = new Matrix(List.transpose(elements))    

    private def dotVectors(a: List[Double],   
                       b: List[Double]): Double = {  
      val multipliedElements =   
        List.map2(a,b)(_*_)  
      (0.0 /: multipliedElements)(_+_)  
    }  

    def *(other: Matrix): Matrix = {  
      require(nCols == other.nRows)  
      val t = other.transpose()  
      new Matrix(  
        for (row <- elements) yield {  
          for (otherCol <- t.elements)  
            yield dotVectors(row, otherCol)  
        }  
      )  

    override def toString(): String = {  
        val rowStrings =   
        for (row <- elements)   
          yield row.mkString("[", ", ", "]")  
        rowStrings.mkString("", "\n", "\n")  
    }  
}  

/* Matrix constructor from a bunch of numbers */
object Matrix {  
  def apply(nCols: Int, els: Double*):Matrix = {  
    def splitRowsWorker(  
      inList: List[Double],   
      working: List[List[Double]]):  
    List[List[Double]] =  
      if (inList.isEmpty)  
        working  
      else {  
        val (a, b) = inList.splitAt(nCols)  
        splitRowsWorker(b, working + a)  
      }  
    def splitRows(inList: List[Double]) =  
      splitRowsWorker(inList, List[List[Double]]())  
    val rows: List[List[Double]] =   
      splitRows(els.toList)  
    new Matrix(rows)  
  }  
}  

EDIT I understood that strictly speaking the answer is No: overloading * is not possible without side-effects of defining also a + and others or special tricks. The numeric-prelude package describes it best:

In some cases, the hierarchy is not finely-grained enough: Operations that are often defined independently are lumped together. For instance, in a financial application one might want a type "Dollar", or in a graphics application one might want a type "Vector". It is reasonable to add two Vectors or Dollars, but not, in general, reasonable to multiply them. But the programmer is currently forced to define a method for '(*)' when she defines a method for '(+)'.

Answer 1

It'll be perfectly safe with a smart constructor and stored dimensions. Of course there are no natural implementations for the operations signum and fromIntegral (or maybe a diagonal matrix would be fine for the latter).

module Matrix (Matrix(),matrix,matrixTranspose) where

import Data.List (transpose)

data Matrix a = Matrix {matrixN :: Int, 
                        matrixM :: Int,
                        matrixElems :: [[a]]}
                deriving (Show, Eq)

matrix :: Int -> Int -> [[a]] -> Matrix a
matrix n m vals
  | length vals /= m            = error "Wrong number of rows"
  | any (/=n) $ map length vals = error "Column length mismatch"
  | otherwise = Matrix n m vals

matrixTranspose (Matrix m n vals) = matrix n m (transpose vals)

instance Num a => Num (Matrix a) where

  (+) (Matrix m n vals) (Matrix m' n' vals')
    | m/=m' = error "Row number mismatch"
    | n/=n' = error "Column number mismatch"
    | otherwise = Matrix m n (zipWith (zipWith (+)) vals vals')

  abs (Matrix m n vals) = Matrix m n (map (map abs) vals)

  negate (Matrix m n vals) = Matrix m n (map (map negate) vals)

  (*) (Matrix m n vals) (Matrix n' p vals')
    | n/=n' = error "Matrix dimension mismatch in multiplication"
    | otherwise = let tvals' = transpose vals'
                      dot x y = sum $ zipWith (*) x y
                      result = map (\col -> map (dot col) tvals') vals
                  in Matrix m p result

Test it in ghci:

*Matrix> let a = matrix 3 2 [[1,0,2],[-1,3,1]]
*Matrix> let b = matrix 2 3 [[3,1],[2,1],[1,0]]
*Matrix> a*b
Matrix {matrixN = 3, matrixM = 3, matrixElems = [[5,1],[4,2]]}

Since my Num instance is generic, it even works for complex matrices out of the box:

Prelude Data.Complex Matrix> let c = matrix 2 2 [[0:+1,1:+0],[5:+2,4:+3]]
Prelude Data.Complex Matrix> let a = matrix 2 2 [[0:+1,1:+0],[5:+2,4:+3]]
Prelude Data.Complex Matrix> let b = matrix 2 3 [[3:+0,1],[2,1],[1,0]]
Prelude Data.Complex Matrix> a
Matrix {matrixN = 2, matrixM = 2, matrixElems = [[0.0 :+ 1.0,1.0 :+ 0.0],[5.0 :+ 2.0,4.0 :+ 3.0]]}
Prelude Data.Complex Matrix> b
Matrix {matrixN = 2, matrixM = 3, matrixElems = [[3.0 :+ 0.0,1.0 :+ 0.0],[2.0 :+ 0.0,1.0 :+ 0.0],[1.0 :+ 0.0,0.0 :+ 0.0]]}
Prelude Data.Complex Matrix> a*b
Matrix {matrixN = 2, matrixM = 3, matrixElems = [[2.0 :+ 3.0,1.0 :+ 1.0],[23.0 :+ 12.0,9.0 :+ 5.0]]}

---

EDIT: new material

Oh, you want to just override the (*) function without any Num stuff. That's possible to o but you'll have to remember that the Haskell standard library has reserved (*) for use in the Num class.

module Matrix where

import qualified Prelude as P
import Prelude hiding ((*))
import Data.List (transpose)

class Multiply a where
  (*) :: a -> a -> a

data Matrix a = Matrix {matrixN :: Int, 
                        matrixM :: Int,
                        matrixElems :: [[a]]}
                deriving (Show, Eq)

matrix :: Int -> Int -> [[a]] -> Matrix a
matrix n m vals
  | length vals /= m            = error "Wrong number of rows"
  | any (/=n) $ map length vals = error "Column length mismatch"
  | otherwise = Matrix n m vals

matrixTranspose (Matrix m n vals) = matrix n m (transpose vals)

instance P.Num a => Multiply (Matrix a) where
  (*) (Matrix m n vals) (Matrix n' p vals')
    | n/=n' = error "Matrix dimension mismatch in multiplication"
    | otherwise = let tvals' = transpose vals'
                      dot x y = sum $ zipWith (P.*) x y
                      result = map (\col -> map (dot col) tvals') vals
                  in Matrix m p result


a = matrix 3 2 [[1,2,3],[4,5,6]]
b = a * matrixTranspose

Testing in ghci:

*Matrix> b
Matrix {matrixN = 3, matrixM = 3, matrixElems = [[14,32],[32,77]]}

There. Now if a third module wants to use both the Matrix version of (*) and the Prelude version of (*) it'll have to of course import one or the other qualified. But that's just business as usual.

I could've done all of this without the Multiply type class but this implementation leaves our new shiny (*) open for extension in other modules.

Answer 2

Alright, there's a lot of confusion about what's happening here floating around, and it's not being helped by the fact that the Haskell term "class" does not line up with the OO term "class" in any meaningful way. So let's try to make a careful answer. This answer starts with Haskell's module system.

In Haskell, when you import a module Foo.Bar, it creates a new set of bindings. For each variable x exported by the module Foo.Bar, you get a new name Foo.Bar.x. In addition, you may:

• import qualified or not. If you import qualified, nothing more happens. If you do not, an additional name without the module prefix is defined; in this case, just plain old x is defined.
• change the qualification prefix or not. If you import as Alias, then the name Foo.Bar.x is not defined, but the name Alias.x is.
• hide certain names. If you hide name foo, then neither the plain name foo nor any qualified name (like Foo.Bar.foo or Alias.foo) is defined.

Furthermore, names may be multiply defined. For example, if Foo.Bar and Baz.Quux both export the variable x, and I import both modules without qualification, then the name x refers to both Foo.Bar.x and Baz.Quux.x. If the name x is never used in the resulting module, this clash is ignored; otherwise, a compiler error asks you to provide more qualification.

Finally, if none of your imports mention the module Prelude, the following implicit import is added:

import Prelude

This imports the Prelude without qualification, with no additional prefix, and without hiding any names. So it defines "bare" names and names prefixed by Prelude., and nothing more.

Here ends the bare basics you need to understand about the module system. Now let's discuss the bare basics you need to understand about typeclasses.

A typeclass includes a class name, a list of type variables bound by that class, and a collection of variables with type signatures that refer to the bound variables. Here's an example:

class Foo a where
    foo :: a -> a -> Int

The class name is Foo, the bound type variable is a, and there is only one variable in the collection, namely foo, with type signature a -> a -> Int. This class declares that some types have a binary operation, named foo, which computes an Int. Any type may later (even in another module) be declared to be an instance of this class: this involves defining the binary operation above, where the bound type variable a is substituted with the type you are creating an instance for. As an example, we might implement this for integers by the instance:

instance Foo Int where
    foo a b = (a `mod` 76) * (b + 7)

Here ends the bare basics you need to understand about typeclasses. We may now answer your question. The only reason the question is tricky is because it falls smack dab on the intersection between two name management techniques: modules and typeclasses. Below I discuss what this means for your specific question.

The module Prelude defines a typeclass named Num, which includes in its collection of variables a variable named *. Therefore, we have several options for the name *:

• If the type signature we desire happens to follow the pattern a -> a -> a, for some type a, then we may implement the Num typeclass. We therefore extend the Num class with a new instance; the name Prelude.* and any aliases for this name are extended to work for the new type. For matrices, this would look like, for example,

  instance Num Matrix where
      m * n = {- implementation goes here -}
  

• We may define a different name than *.

  m |*| n = {- implementation goes here -}
  

• We may define the name *. Whether this name is defined as part of a new type class or not is immaterial. If we do nothing else, there will then be at least two definitions of *, namely, the one in the current module and the one implicitly imported from the Prelude. We have a variety of ways of dealing with this. The simplest is to explicitly import the Prelude, and ask for the name * not to be defined:

  import Prelude hiding ((*))
  

  You might alternately choose to leave the implicit import of Prelude, and use a qualified * everywhere you use it. Other solutions are also possible.

The main point I want you to take away from this is: the name * is in no way special. It is just a name defined by the Prelude, and all of the tools we have available for namespace control are available.

Answer 3

You can implement * as matrix multiplication by defining an instance of Num class for Matrix. But the code won't be type-safe: * (and other arithmetic operations) on matrices as you define them is not total, because of size mismatch or in case of '/' non-existence of inverse matrices.

As for 'the hierarchy is not defined precisely' - there is also Monoid type class, exactly for the cases when only one operation is defined.

There are too many things to be 'added', sometimes in rather exotic ways (think of permutation groups). Haskell designers designed to reserve arithmetical operations for different representations of numbers, and use other names for more exotic cases.