Matrix Operations

June 22, 2010

All of these operations rely heavily on loops using the for syntax given in the Standard Prelude. Here is addition:

(define (matrix-add a b)
  (let ((ar (matrix-rows a)) (ac (matrix-cols a))
        (br (matrix-rows b)) (bc (matrix-cols b)))
    (if (or (not (= ar br)) (not (= ac bc)))
        (error 'matrix-add "incompatible matrices")
        (let ((c (make-matrix ar ac)))
          (for (i ar)
            (for (j ac)
              (matrix-set! c i j
                (+ (matrix-ref a i j)
                   (matrix-ref b i j)))))
          c))))

Multiplying a matrix by a scalar is just as simple. Note that the scalar, as well as the matrix elements, can be any numeric type:

(define (matrix-scalar-multiply n a)
  (let* ((ar (matrix-rows a))
         (ac (matrix-cols a))
         (c (make-matrix ar ac)))
    (for (i ar)
      (for (j ac)
        (matrix-set! c i j
          (* n (matrix-ref a i j)))))
    c))

Multiplying two matrices involves a third for loop:

(define (matrix-multiply a b)
  (let ((ar (matrix-rows a)) (ac (matrix-cols a))
        (br (matrix-rows b)) (bc (matrix-cols b)))
    (if (not (= ac br))
        (error 'matrix-multiply "incompatible matrices")
        (let ((c (make-matrix ar bc 0)))
          (for (i ar)
            (for (j bc)
              (for (k ac)
                (matrix-set! c i j
                  (+ (matrix-ref c i j)
                     (* (matrix-ref a i k)
                        (matrix-ref b k j)))))))
          c))))

Transposing a matrix involves no arithmetic at all, just subscript manipulations:

(define (matrix-transpose a)
  (let* ((ar (matrix-rows a))
         (ac (matrix-cols a))
         (c (make-matrix ac ar)))
    (for (i ar)
      (for (j ac)
        (matrix-set! c j i
          (matrix-ref a i j))))
    c))

Here are some examples:

> (define a #(#(1)
              #(2)
              #(3)))
> (define b #(#(1 2 3)
              #(4 5 6)))
> (define c #(#(2 3 4)
              #(3 4 5)))
> (define d #(#(1 2 3 4)
              #(2 3 4 5)
              #(3 4 5 6)))
> (matrix-add b c)
#2(#3(3 5 7)
   #3(7 9 11))
> (matrix-scalar-multiply 2 b)
#2(#3(2 4 6)
   #3(8 10 12))
> (matrix-multiply b d)
#2(#4(14 20 26 32)
   #4(32 47 62 77))
> (matrix-transpose b)
#3(#2(1 4)
   #2(2 5)
   #2(3 6))

We used the matrix operations and for operator from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/5nzWMChE.

About these ads

Pages: 1 2

12 Responses to “Matrix Operations”

  1. [...] Praxis – Matrix Operations By Remco Niemeijer In today’s Programming Praxis exercise our task is to implement four common matrix operations. The provided [...]

  2. Remco Niemeijer said

    My Haskell solution (see http://bonsaicode.wordpress.com/2010/06/22/programming-praxis-matrix-operations/ for a version with comments):

    add :: Num a => [[a]] -> [[a]] -> [[a]]
    add = zipWith $ zipWith (+)
    
    scale :: Num a => a -> [[a]] -> [[a]]
    scale = map . map . (*)
    
    transpose :: [[a]] -> [[a]]
    transpose [] = []
    transpose xs = foldr (zipWith (:)) (repeat []) xs
    
    mult :: Num a => [[a]] -> [[a]] -> [[a]]
    mult a b = [map (sum . zipWith (*) r) $ transpose b | r <- a]
    
  3. Here’s my Java solution.

    @Remco,

    Your Haskell solution is gorgeous (its simplicity astounds me). I’ll hopefully be able to churn out solutions like that in Haskell soon! =)

  4. Remco Niemeijer said

    @Elvis Montero:

    Thanks. Best of luck with learning Haskell. I can highly recommend it :)

  5. Paul Carleton said

    @Remco,

    Your post has inspired me to learn Haskell!

  6. Geir S said

    ;;; A scheme solution. Nothing like the beauty of the haskell one though ;-)

    ;; Q1
    (define (add-matrix m1 m2)
    (cond ((or (null? m1) (null? m2)) ‘())
    (else (cons (sum-list (car m1) (car m2))
    (add-matrix (cdr m1) (cdr m2))))))

    ;; Q1 helper
    (define (sum-list l1 l2)
    (cond ((or (null? l1) (null? l2)) ‘())
    (else (cons (+ (car l1) (car l2))
    (sum-list (cdr l1) (cdr l2))))))

    ;; Q2
    (define (scale s m)
    (letrec ((scale* (lambda (m)
    (cond ((null? m) ‘())
    (else (cons (mul-list s (car m))
    (scale* (cdr m))))))))
    (scale* m)))

    ;; Q2 helper
    (define (mul-list s l)
    (letrec ((mul-list* (lambda (l)
    (cond ((null? l) ‘())
    (else (cons (* s (car l))
    (mul-list* (cdr l))))))))
    (mul-list* l)))

    ;; Q3
    (define (X A B)
    (cond ((null? A) ‘())
    (else (cons (psum (car A) B)
    (X (cdr A) B)))))

    ;; Q3 helper – genereate list of cross-sums
    (define (psum l m)
    (letrec ((psum* (lambda ™
    (cond ((null? tm) ‘())
    (else (cons (xs l (car tm))
    (psum* (cdr tm))))))))
    (psum* (transpose m))))

    ;; Q3 helper – multiply each element and sum result
    (define (xs l1 l2)
    (cond ((or (null? l1) (null? l2)) 0)
    (else (+ (* (car l1) (car l2)) (xs (cdr l1) (cdr l2))))))

    ;; Q4
    (define (transpose m)
    (cond ((null? (car m)) ‘())
    (else (cons (list-of-heads m)
    (transpose (list-of-tails m))))))

    ;; Q4 helper
    (define (list-of-heads m)
    (cond ((null? m) ‘())
    (else (cons (car (car m))
    (list-of-heads (cdr m))))))

    ;; Q4 helper
    (define (list-of-tails m)
    (cond ((null? m) ‘())
    (else (cons (cdr (car m))
    (list-of-tails (cdr m))))))

  7. Geir S said

    bah.. indents lost :-(

    one with proper indents on pocoo….

    http://paste.pocoo.org/show/nZhGWlzvWmDhV9pRV5im/

  8. erislover said

    Geir S: list-of-heads is just (map car m)

  9. erislover said

    Interestingly I realized after writing it that I duplicated the Haskell version above in Scheme, though I put scalar multiplication into matrix multiplication. I have such a hard time reading Haskell that I couldn’t notice it before! No error checking. A matrix is a list of lists.

    (define transpose
      (λ(l)
        (if (empty? (first l))
            empty
            (cons (map first l)
                  (transpose (map rest l))))))
    
    (define add
      (λ(one two)
        (map (λ(x y)
               (map + x y))
             one two)))
    
    (define times
      (λ(one two)
        (let ((sum (λ(list)
                     (foldl + 0 list))))
          (if (number? one)
              (map (λ(x) (map (λ(y) (* one y)) x)) two)
              (let ((two (transpose two)))
                (map (λ(x) (map (λ(y) (sum (map * x y))) two)) one))))))
    
  10. Eric Pierce said

    My Haskell solution:

    -- Matrix.hs
    --
    -- Defines a Matrix data type and operations upon Matrices
    
    data Matrix a = Matrix [[a]]
         deriving (Show, Eq)
    
    -- Add two matrices together
    --
    --   > addMatrices (Matrix [[1,2,3], [4,5,6]]) (Matrix [[2,3,4], [3,4,5]])
    --   Matrix [[3,5,7], [7,9,11]]
    --
    addMatrices :: (Num a) => Matrix a -> Matrix a -> Matrix a
    addMatrices (Matrix m1) (Matrix m2) = Matrix $ zipWith (zipWith (+)) m1 m2 
    
    -- Multiply a matrix by a scalar
    --
    --   > multScalar 2 (Matrix [[1,2,3], [4,5,6]])
    --   Matrix [[2,4,6], [8,10,12]]
    --
    multScalar :: (Num a) => a -> Matrix a -> Matrix a
    multScalar x (Matrix m) = Matrix $ map (map (x*)) m
    
    -- Multiply two matrices together
    --
    --   > multMatrices (Matrix [[1,2,3], [4,5,6]]) (Matrix [[1,2,3,4], [2,3,4,5], [3,4,5,6]])
    --   Matrix [[14,20,26,32],[32,47,62,77]]
    --
    multMatrices :: (Num a) => Matrix a -> Matrix a -> Matrix a
    multMatrices (Matrix m1) (Matrix m2) = Matrix $ [ map (multRow r) m2t | r <- m1 ]
        where (Matrix m2t) = transposeMatrix (Matrix m2)
              multRow r1 r2 = sum $ zipWith (*) r1 r2 
    
    -- Transpose a matrix
    --
    --   > transposeMatrix (Matrix [[1,2,3], [4,5,6]])
    --   Matrix [[1,4], [2,6], [3,6]]
    --
    transposeMatrix :: Matrix a -> Matrix a
    transposeMatrix (Matrix m) = Matrix (zipList m) 
    
    -- Zip together a list of lists. The result is truncated to the length of the
    -- shortest list. This is like the builtin zip function, except it can zip
    -- together an arbitrary number of lists.
    zipList :: [[a]] -> [[a]]
    zipList lists
        -- Base cases: there are no lists, or any sub-list is empty
        | length lists == 0          = [] 
        | any ((==0) . length) lists = [] 
        -- Take the head from each sub-list and recurse with the tails
        | otherwise                  = map head lists : zipList (map tail lists)
    

    Not nearly as sexy as @Remco’s, partly because I used a Matrix data type that complicated things somewhat unnecessarily. I wrote zipList from scratch before remembering there was a ‘transpose’ function in Data.List, but overall this was a good learning experience. I remember back when I was into C++, and I wrote some matrix classes–the header declarations alone were longer than this :-)

  11. My Python solution with convenient operators overloading, though not very effective:

    class Matrix:
    	def __init__(self, rows):
    		if not isinstance(rows, list) or \
    		not isinstance(rows[0], list) or \
    		len(rows) == 0 or len(rows[0]) == 0:
    			raise ValueError
    		self.rows, self.cols = len(rows), len(rows[0])
    		for row in rows:
    			if len(row) != self.cols:
    				raise ValueError
    		self.m = list(map(list, rows))
    		
    	def __repr__(self):
    		return '\n'.join(map(str, self.m))
    	
    	def __add__(self, other):
    		data = []
    		for row1, row2 in zip(self.m, other.m):
    			data.append(list(map(lambda a, b: a + b, row1, row2)))
    		return Matrix(data)
    	
    	def __mul__(self, other):
    		if isinstance(self, Matrix) and isinstance(other, Matrix):
    			return self.matrix_mult(other)
    		else:
    			return self.scalar_mult(other)
    
    	def __rmul__(self, other):
    		if isinstance(self, Matrix) and isinstance(other, Matrix):
    			return self.matrix_mult(other)
    		else:
    			return self.scalar_mult(other)
    	
    	def matrix_mult(self, other):
    		if self.cols != other.rows:	
    			raise ValueError
    		data = []
    		for i in range(self.rows):
    			row = []
    			for j in range(other.cols):
    				row.append(sum(map(lambda a: a[0]*a[1], zip(self.m[i], map(lambda l: l[j], other.m)))))
    			data.append(row)
    		return Matrix(data)
    	
    	def scalar_mult(self, scalar):
    		data = []
    		for row in self.m:
    			data.append(list(map(lambda x: x * scalar, row)))
    		return Matrix(data)
    	
    	def transpose(self):
    		data = []
    		for i in range(self.cols):
    			data.append(list(map(lambda x: x[i], self.m)))
    		return Matrix(data)
    
  12. Will said

    Scheme solution:

    (define (add mat1 mat2) (map (lambda (x y) (map + x y)) mat1 mat2))
    
    (define (scale n mat) (map (lambda (y) (map (lambda (x) (* n x)) y)) mat))
    
    (define (transpose mat) (apply map list mat))
    
    (define (mult mat1 mat2)
      (map
       (lambda (row) (apply map (lambda column (apply + (map * row column))) mat2))
       mat1))
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 576 other followers

%d bloggers like this: