## Sunday, April 22, 2012

### SICP 2.56 - 2.58: Symbolic Differentiation

From SICP section 2.3.2  Example: Symbolic Differentiation

This section of SICP presents an application of symbol manipulation in the form of the following procedure that performs symbolic differentiation of algebraic expressions.

(define (deriv exp var)
(cond ((number? exp) 0)
((variable? exp)
(if (same-variable? exp var) 1 0))
((sum? exp)
(make-sum (deriv (addend exp) var)
(deriv (augend exp) var)))
((product? exp)
(make-sum
(make-product (multiplier exp)
(deriv (multiplicand exp) var))
(make-product (deriv (multiplier exp) var)
(multiplicand exp))))
(else
(error "unknown expression type -- DERIV" exp))))


The following supporting procedures are required by deriv. Each one is explained in the text.

(define (variable? x) (symbol? x))

(define (same-variable? v1 v2)
(and (variable? v1) (variable? v2) (eq? v1 v2)))

(define (make-sum a1 a2) (list '+ a1 a2))
(define (make-product m1 m2) (list '* m1 m2))

(define (sum? x)
(and (pair? x) (eq? (car x) '+)))

(define (addend s) (cadr s))
(define (augend s) (caddr s))

(define (product? x)
(and (pair? x) (eq? (car x) '*)))

(define (multiplier p) (cadr p))
(define (multiplicand p) (caddr p))


With these definitions, we can perform differentiation on algebraic expressions in Scheme's familiar prefix notation.

> (deriv '(+ x 3) 'x)
(+ 1 0)
> (deriv '(* x y) 'x)
(+ (* x 0) (* 1 y))
> (deriv '(* (* x y) (+ x 3)) 'x)
(+ (* (* x y) (+ 1 0)) (* (+ (* x 0) (* 1 y)) (+ x 3)))


This procedure returns results that are correct, but the expressions are not in their simplest form. For example, instead of (+ 1 0) we'd like the procedure to simplify the expression to 1. Similarly, the expression (+ (* x 0) (* 1 y)) can be simplified to y. We only need to change the constructors make-sum and make-product to simplify the resulting expressions. We can change make-sum so that if both its arguments are numbers, make-sum will return their sum. Also, if either of the arguments is 0, then make-sum will return the other summand (whether it's a number or a symbol).

(define (make-sum a1 a2)
(cond ((=number? a1 0) a2)
((=number? a2 0) a1)
((and (number? a1) (number? a2)) (+ a1 a2))
(else (list '+ a1 a2))))


Similarly, make-product should return 0 if either of its arguments is 0, and if either of its arguments is 1 it should return the other argument.

(define (make-product m1 m2)
(cond ((or (=number? m1 0) (=number? m2 0)) 0)
((=number? m1 1) m2)
((=number? m2 1) m1)
((and (number? m1) (number? m2)) (* m1 m2))
(else (list '* m1 m2))))


Both of these new implementations make use of the =number? procedure to check whether an expression is equal to an expected value.

(define (=number? exp num)
(and (number? exp) (= exp num)))


Here's how the improved procedures work with the same three examples.

> (deriv '(+ x 3) 'x)
1
> (deriv '(* x y) 'x)
y
> (deriv '(* (* x y) (+ x 3)) 'x)
(+ (* x y) (* y (+ x 3)))


The first two resulting expressions are now in simplest form. The third is greatly improved, but could still be better.

Exercise 2.56 asks us to extend the deriv procedure to handle more kinds of expressions. We're to implement the differentiation rule

$\frac{d(u^n)}{dx} = nu^{(n-1)} (\frac{du}{dx})$

by adding a new clause to the deriv program and defining appropriate procedures exponentiation?, base, exponent, and make-exponentiation. (We'll use the symbol ** to denote exponentiation.) We'll also build in the rules that anything raised to the power 0 is 1 and any expression raised to the power 1 is the expression itself. We can base the implementations of exponentiation?, base, and exponent on the corresponding procedures used for sums and products.

(define (exponentiation? x)
(and (pair? x) (eq? (car x) '**)))

(define (base e) (cadr e))
(define (exponent e) (caddr e))

Similarly, our implementation of make-exponentiation will be based on the updated versions of make-sum and make-product, since we're going to build in two rules of exponentiation that will simplify resulting expressions.

(define (make-exponentiation base exponent)
(cond ((=number? exponent 0) 1)
((=number? exponent 1) base)
((and (number? base) (number? exponent)) (expt base exponent))
(else (list '** base exponent))))


The last step before testing is to extend the deriv procedure itself to recognize and handle exponentiation.

(define (deriv exp var)
(cond ((number? exp) 0)
((variable? exp)
(if (same-variable? exp var) 1 0))
((sum? exp)
(make-sum (deriv (addend exp) var)
(deriv (augend exp) var)))
((product? exp)
(make-sum
(make-product (multiplier exp)
(deriv (multiplicand exp) var))
(make-product (deriv (multiplier exp) var)
(multiplicand exp))))
((exponentiation? exp)
(make-product
(make-product (exponent exp)
(make-exponentiation
(base exp)
(make-sum (exponent exp) -1)))
(deriv (base exp) var)))
(else
(error "unknown expression type -- DERIV" exp))))


We can test the extended implementation of deriv using an example from the beginning of the current section of the text. "For example, if the arguments to the procedure are $ax^2 + bx + c$ and x, the procedure should return $2ax + b$." I'll translate the expression to prefix notation in multiple steps so we can see the effect that each term has on the result.

> (deriv '(** x 2) 'x)
(* 2 x)
> (deriv '(* a (** x 2)) 'x)
(* a (* 2 x))
> (deriv '(+ (* a (** x 2)) (* b x)) 'x)
(+ (* a (* 2 x)) b)
> (deriv '(+ (+ (* a (** x 2)) (* b x)) c) 'x)
(+ (* a (* 2 x)) b)


Exercise 2.57 asks us to extend the deriv program to handle sums and products with arbitrary numbers of (two or more) terms. The third example above could be expressed as

(deriv '(* x y (+ x 3)) 'x)

We'll do this by changing only the representation for sums and products, without changing the deriv procedure itself. (For example, the addend of a sum would be the first term, and the augend would be the sum of the rest of the terms.) We can change the representation for sums and products by changing only two procedures, augend and multiplicand. The original implementations of these procedures simply returned the second argument in a sum or product, respectively. Instead of returning a single value or symbol, we'll need to modify these procedures so they can return an expression.

(define (augend s)
(if (null? (cdddr s))
(caddr s)
(cons '+ (cddr s))))


The new augend procedure first tests to see if the expression passed in has only two arguments by checking to see if its third argument is null. If the expression has only two arguments, the second one is returned, just as it would have been in the original implementation. Otherwise, if the expression has more than two arguments, all of the arguments following the first one are returned in a new sum expresssion. The new implementation of multiplicand is pretty much the same.

(define (multiplicand p)
(if (null? (cdddr p))
(caddr p)
(cons '* (cddr p))))


We can use the example from the text to test that the deriv procedure works correctly with the new representations of sums and products.

> (deriv '(* x y (+ x 3)) 'x)
(+ (* x y) (* y (+ x 3)))


Exercise 2.58 asks us to modify the differentiation program so that it works with ordinary mathematical notation (where + and * are infix rather than prefix operators). Since the differentiation program is defined in terms of abstract data, we can modify it to work with different representations of expressions solely by changing the predicates, selectors, and constructors that define the representation of the algebraic expressions on which the differentiator is to operate. We're to do this in two steps, a fully parenthesized version and a version that drops unnecessary parentheses:

a. Show how to do this in order to differentiate algebraic expressions presented in infix form, such as (x + (3 * (x + (y + 2)))). To simplify the task, assume that + and * always take two arguments and that expressions are fully parenthesized.

We can solve the first part of the problem simply by changing the procedures that define how a sum is represented.  Instead of the + symbol appearing first in an expression, it will now appear second.

(define (make-sum a1 a2)
(cond ((=number? a1 0) a2)
((=number? a2 0) a1)
((and (number? a1) (number? a2)) (+ a1 a2))
(else (list a1 '+ a2))))

(define (sum? x)
(and (pair? x) (eq? (cadr x) '+)))

(define (addend s) (car s))
(define (augend s) (caddr s))


The definitions for products are equivalent.

(define (make-product m1 m2)
(cond ((or (=number? m1 0) (=number? m2 0)) 0)
((=number? m1 1) m2)
((=number? m2 1) m1)
((and (number? m1) (number? m2)) (* m1 m2))
(else (list m1 '* m2))))

(define (product? x)
(and (pair? x) (eq? (cadr x) '*)))

(define (multiplier p) (car p))
(define (multiplicand p) (caddr p))


We can test with a few simple examples before moving to the more complicated example given in the text.

> (deriv '(x + 3) 'x)
1
> (deriv '(x * y) 'x)
y
>  (deriv '(x + (3 * (x + (y + 2)))) 'x)
4


b. The problem becomes substantially harder if we allow standard algebraic notation, such as (x + 3 * (x + y + 2)), which drops unnecessary parentheses and assumes that multiplication is done before addition. Can you design appropriate predicates, selectors, and constructors for this notation such that our derivative program still works?

Only a few additional changes are necessary in order to correctly interpret expressions where unnecessary parentheses are excluded. In part a above we defined both the augend and multiplicand of an expression using caddr. We were able to do this because we knew all expressions would be two parameters separated by an operator, and that the expression would be contained in parentheses.

Since this is no longer the case, we now want the augend and multiplicand procedures to return an entire subexpression, so we'd like to use cddr. The only problem with this is that cddr always returns a list, so it will fail in those cases where the augend or multiplicand is a single value or symbol, as would be the case in fully-parenthesized expressions. We can get around this limitation by introducing a procedure that will simplify sub-expressions of this form by returning a single value or symbol if that's all there is to an expression.

(define (simplify exp)
(if (null? (cdr exp)) (car exp) exp))


Now we can define augend and multiplicand in terms of simplify and cddr.

(define (augend s) (simplify (cddr s)))
(define (multiplicand p) (simplify (cddr p)))


We can test a few different expressions, starting with the example given.

> (deriv '(x + 3 * (x + y + 2)) 'x)
4
> (deriv '(x + 3 * (x + y + 2)) 'y)
3


Finally, we should also test a few simpler cases to make sure our changes didn't break anything.

> (deriv '(x + 3) 'x)
1
> (deriv '(x * y) 'x)
y
> (deriv '(x * y * (x + 3)) 'x)
((x * y) + (y * (x + 3)))


Related:

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

#### 2 comments:

Alexey Grigorev said...

Hello!

Try running (deriv '(x * 3 + 5 * (x + y + 2)) 'x) for the last task

Anonymous said...

Hey Bill, your solution for 2.58.b doesn't work correctly.
For the example given by Alexey above, the derivative of (x * 3 + 5 * (x + y + 2)) with respect to x is 8 but your code gives an incorrect result.
Here is a solution for 2.58.b that works correctly.
I got this from several solutions at http://community.schemewiki.org/?sicp-solutions

(define (operation expr)
(cond ((memq '+ expr) '+)
((memq '* expr) '*)
((memq '** expr) '**)))

(define (prefix item x)
(define (iter x result)
(if (eq? (car x) item)
(reverse result)
(iter (cdr x) (cons (car x) result))))
(iter x '()))

(define (sum? x) (eq? '+ (operation x)))

(define (addend s)
(let ((a (prefix '+ s)))
(if (= (length a) 1) (car a) a)))

(define (augend s)
(let ((a (cdr (memq '+ s))))
(if (= (length a) 1) (car a) a)))

(define (product? x) (eq? '* (operation x)))

(define (multiplier p)
(let ((m (prefix '* p)))
(if (= (length m) 1) (car m) m)))

(define (multiplicand p)
(let ((m (cdr (memq '* p))))
(if (= (length m) 1) (car m) m)))

(define (exponentiation? x) (eq? '** (operation x)))

(define (base e)
(let ((p (prefix '** e)))
(if (= (length p) 1) (car p) p)))

(define (exponent e)
(let ((p (cdr (memq '** e))))
(if (= (length p) 1) (car p) p)))