Sunday, April 24, 2011

SICP 2.36 & 2.37: Matrix Algebra

From SICP section 2.2.3 Sequences as Conventional Interfaces

Exercise 2.36 asks us to complete the accumulate-n procedure, which is similar to accumulate except that it takes as its third argument a sequence of sequences that are assumed to all have the same length. It applies the accumulation procedure to all the first elements of the sub-sequences, all the second elements, third elements, and so on and returns a sequence of the results. For example, given the sequence s containing the following values, ((1 2 3) (4 5 6) (7 8 9) (10 11 12)), the value of (accumulate-n + 0 s) should return the sequence (22 26 30).

Our task is to fill in the missing expressions in the following definition:
(define (accumulate-n op init seqs)
(if (null? (car seqs))
null
(cons (accumulate op init <??>)
(accumulate-n op init <??>))))

If you've been doing all the exercises up to this point, then this procedure's overall structure of recursively building up a sequence with cons will be familiar.

In exercise 2.35 we used map and a simple lambda expression to flatten out a tree so the result could be used as the third parameter to accumulate. We can do something very similar here to fill in the missing expressions above. We'll still use map to loop over each of the sub-sequences, but in this case we need map to return a sequence that contains the first element from each sub-sequence (in the first missing expression) and the remaining elements of each sub-sequence (in the second missing expression). We already know two procedures that meet those needs, car and cdr.
; accumulate from sicp section 2.2.3
(define (accumulate op initial sequence)
(if (null? sequence)
initial
(op (car sequence)
(accumulate op initial (cdr sequence)))))

(define (accumulate-n op init seqs)
(if (null? (car seqs))
null
(cons (accumulate op init (map car seqs))
(accumulate-n op init (map cdr seqs)))))

The first call to map returns a sequence containing the first element of each of the original sub-sequences. The second call to map returns the sequence of sub-sequences with each of their first elements removed. This continues recursively until no more elements remain. We can test the solution with the values given in the exercise.
> (define s (list (list 1 2 3) (list 4 5 6) (list 7 8 9) (list 10 11 12)))
> s
((1 2 3) (4 5 6) (7 8 9) (10 11 12))
> (accumulate-n + 0 s)
(22 26 30)


Exercise 2.37 introduces a representation for vectors and matrices. A vector can be represented as a sequence of numbers, and a matrix as a sequence of vectors. For example, the matrix


can be represented as the sequence ((1 2 3 4) (4 5 6 6) (6 7 8 9)).

We can use this representation to express the basic operations of matrix algebra in terms of sequence operations. These basic operations are:

(dot-product v w) - returns the sum Σi(viwi)

(matrix-*-vector m v) - returns the vector t, where ti = Σj(mijvj)

(matrix-*-matrix m n) - returns the matrix p, where pij = Σk(miknkj)

(transpose mat) - returns the matrix n, where nij = mji


Dot product

The dot product is defined for us.
(define (dot-product v w)
(accumulate + 0 (map * v w)))

The dot-product procedure takes two vectors (that are assumed to be equal length) and returns a single number obtained by multiplying corresponding elements and then summing those products.

Note that the use of map in this procedure takes more arguments than we've seen up to this point. That's because map takes a procedure of n arguments, together with n lists, and applies the procedure to all the first elements of the lists, all the second elements of the lists, and so on, returning a list of the results. Up to this point, we've been using map where n = 1.

Our task is to fill in the missing expressions in the given procedures for computing the remaining matrix operations. Let's take them one at a time.


Multiply a matrix by a vector
(define (matrix-*-vector m v)
(map <??> m))

When multiplying a matrix by a vector, each row of the matrix is multiplied by the vector (using the dot product procedure described above). The result is a vector containing each dot product. Since map will pass each row of the matrix to whatever function we provide, we can just define a lambda expression that uses the dot-product procedure already defined.
(define (matrix-*-vector m v)
(map (lambda (row) (dot-product row v)) m))

Transpose a matrix
(define (transpose mat)
(accumulate-n <??> <??> mat))

There are several ways of looking at matrix transposition, but the one that makes the most sense given the code above is to write the columns of the input matrix as the rows of the output matrix. Remember that accumulate-n will combine the columns with whatever function we give it, and this part of the exercise is a snap.
(define (transpose mat)
(accumulate-n cons null mat))

We can test this procedure with the matrix s defined in the exercise above.
> (define s (list (list 1 2 3) (list 4 5 6) (list 7 8 9) (list 10 11 12)))
> (transpose s)
((1 4 7 10) (2 5 8 11) (3 6 9 12))


Multiply two matrices
(define (matrix-*-matrix m n)
(let ((cols (transpose n)))
(map <??> m)))

When multiplying two matrices m and n, the resulting matrix will have the same number of rows as m and the same number of columns as n. Each element of the result matrix can be found by taking the dot product of each row of m and each column of n.

Matrix Multiplication
When viewed this way, it's easy to see that each row the resulting matrix is the product of a row in the first input matrix (m) and the entire second input matrix (n). Using this knowledge, we can define multiplication of two matrices in terms of matrix-*-vector.
(define (matrix-*-matrix m n)
(let ((cols (transpose n)))
(map (lambda (row) (matrix-*-vector cols row)) m)))


Testing

We have to remember a few rules of matrix algebra that aren't enforced in the procedures we've defined.
  • For dot-product, two vectors can only be multiplied if they're the same length.
  • For matrix-*-vector, the vector must be the same length as the number of columns in the matrix.
  • For matrix-*-matrix, two matrices A and B can be multiplied only if the number of columns of A matches the number of rows of B.

> (define v (list 1 3 -5))
> (define w (list 4 -2 -1))
> (dot-product v w)
3
> (define m (list (list 1 2 3) (list 4 5 6)))
> (define v (list 1 2 3))
> (matrix-*-vector m v)
(14 32)
> (define a (list (list 14 9 3) (list 2 11 15) (list 0 12 17) (list 5 2 3)))
> (define b (list (list 12 25) (list 9 10) (list 8 5)))
> (matrix-*-matrix a b)
((273 455) (243 235) (244 205) (102 160))

Note that the values above are the same as those used as as examples on pages I linked to explaining each operation. The results above are correct, but you can do a Google search for "matrix multiplication calculator" to find several online calculators if you want to test with other values.


Related:

For links to all of the SICP lecture notes and exercises that I've done so far, see The SICP Challenge.

3 comments:

chege said...

using map in 2.2.3 is all sorts of beautiful. I'v been pulling hair out for hours thinking how to use car and cdr to get the corresponding n-th elements of all lists to be passed to accumulate. I think i am still thinking imperatively, since i even first did it in python with for loops, just to make sure it works.

Anonymous said...

@chege You are not alone:) I think it is the biggest problem in embracing functional paradigm for us - to get rid of imperative thinking. I have the same problem.

Unknown said...

Now I am starting reading SICP and so far got the exerciese 2.37, still seems happy with DrRacket with sicp plugin(https://docs.racket-lang.org/sicp-manual/) but it's still not 100% scheme. For example, in this exercise I am using:
(define (transpose mat)
(accumulate-n cons
(map (lambda (x) (list)) (car mat))
mat))

instead using null directly, otherwise it will complains "null: unbound identifier in module in: null"

What are you guys using for SICP exercise, I am using Linux myself, but Windows in work, and was reasonable happy with emacs with Mit Scheme I think? But end up using DrRacket as it seems slightly eaiser to setup for Linux and (expecially) Windows.

Thanks
Paul