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.

Posted by programmingpraxis

Filed in Exercises

12 Comments »

12 Responses to “Matrix Operations”

Programming Praxis – Matrix Operations « Bonsai Code said
June 22, 2010 at 11:13 AM
[…] Praxis – Matrix Operations By Remco Niemeijer In today’s Programming Praxis exercise our task is to implement four common matrix operations. The provided […]

Remco Niemeijer said

June 22, 2010 at 11:13 AM

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]

Elvis Montero said
June 26, 2010 at 9:45 PM
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! =)
Remco Niemeijer said
June 27, 2010 at 10:06 AM
@Elvis Montero:

Thanks. Best of luck with learning Haskell. I can highly recommend it :)
Paul Carleton said
June 27, 2010 at 11:07 PM
@Remco,

Your post has inspired me to learn Haskell!
Geir S said
June 28, 2010 at 9:36 AM
;;; 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))))))
Geir S said
June 28, 2010 at 9:37 AM
bah.. indents lost :-(

one with proper indents on pocoo….

http://paste.pocoo.org/show/nZhGWlzvWmDhV9pRV5im/
erislover said
July 6, 2010 at 4:43 AM
Geir S: list-of-heads is just (map car m)

erislover said

July 6, 2010 at 7:01 PM

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

Eric Pierce said

July 21, 2010 at 8:33 PM

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

Dmitriy Borodiy said

December 6, 2011 at 2:42 PM

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)

Will said

July 22, 2013 at 9:37 PM

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

Programming Praxis