Derivative, dual numbers and Guile

We learn in high school the differentiation rules which reads,

Nothing fancy!

Although old, dual numbers are not well-know, from what I have seen. Let define a dual number \(x\) with two real numbers \(v\) and \(d\), such that,

and similarly \(y=(u, p)\) another dual number. The addition, multiplication, subtraction and division,

Wait! It looks really similar, well even, the same identical with the differentiation rules. Indeed, let denote this number as \((\texttt{value}, \texttt{deriv})\) where,

1. value behaves as real numbers,
2. deriv behaves as derivatives.

If we think as a C programmer, we would like to use some struct to represent these dual numbers; something like,

typedef struct {
    char name[MAXCHAR];
    float value;
    float deriv;
} dual;

In Guile, structures are named record. And a record is the mechanism to define new disjoint type. Guile provides a syntax for defining new record but it can be extended adding some facilities when loading the module SRFI-9. Yeah, the names based on number is impractical–the best is to rely on this page of the Guile manual.

The struct reads,

(use-modules (srfi srfi-9))

(define-record-type <dual>              ;1 name
  (dual name value deriv)               ;2 constructor
  dual?                                 ;3 predicate
  (name  dual-name)                     ;4 field and accessor
  (value dual-value)                    ;5 idem
  (deriv  dual-deriv))                  ;6 idem

where, some explanations,

1. <dual> is the name of the record; it can be anything and <> is a convention,
2. (dual name value deriv) is the constructor; to create a new dual number, it reads, e.g., (define z (dual "whatever" 1 2)),
3. dual? is a predicate which tests if it is or not the type of <dual>; for instance, (dual? z) return #t (true) whereas (dual? 42) returns #f (false)
4. (name dual-name) where name stores the value provided in the constructor and dual-name is the accessor of that field name,
5. Idem for value,
6. Idem for deriv.

scheme@(guile-user)> (define z (dual "whatever" 1 2))
scheme@(guile-user)> z
$3 = #<<dual> name: "whatever" value: 1 deriv: 2>
scheme@(guile-user)> (dual? z)
$4 = #t
scheme@(guile-user)> (dual? 42)
$5 = #f
scheme@(guile-user)> (dual-name z)
$6 = "whatever"
scheme@(guile-user)> (dual-value z)
$7 = 1

Now, we would like to show a dual number, for instance, get this output,

scheme@(guile-user)> (dual-show z)
"whatever"
Value: 1.0E+0
Deriv: 2.0E+0
$8 = #<<dual> name: "whatever" value: 1 deriv: 2>

For one part, it is simply the use of the accessors. For the other part, it is format – similar to printf from C, though escape sequences are marked with ~ instead of %. The procedure dual-show reads,

(define (dual-show z)
  (format #t "~s~%Value: ~e~%Deriv: ~e~%"
          (dual-name z)
          (dual-value z)
          (dual-deriv z))
  z)

We return the argument for returning something. That's arbitrary at this point.

Previously, we used the dual number \((a, 1.)\), so we would like a procedure to build it. From what we know, we would write,

(define (dual-with-one name val)
  (dual name val 1.))

but instead, let introduce the let binding which allows to declare local variable bindings,

(define (variable name val)
  (let ((local 1.0))
    (dual name val local)))

Here, local is a local variable where the value is set to 1.. Note, conditional as if could be there, for instance,

(let  ((local (if name
                  name
                  "c")))
  ...)

which means, if name is set, then use this value for the variable local, else use the string "c". Hum, we would like to define a new procedure but with an optional argument. It is provided by the variant define*, for instance,

(define* (constant val #:optional name)
  (let  ((local (if name
                    name
                    "c")))
    (dual local val 0)))

where #:optional name means it is… optional. Calling (constant 2) set the default value #f (false) for name. Let as an exercise to set another default value than #f, for example 42. A quick remark: atoms with #: as #:optinonal are named keyword.

Another quick remark: append a star often means it is a variant; define* is an example, it is a variant of define. For instance, let* is also a variant of let: it allows to create local variables which depend on previous declared ones,

(let* ((one 1)
       (two (+ 1 one))
       (three (1+ two)))
  (* 2 three))

Note that 1+ is not a typo.

What is still missing is how to manipulate strings. For instance, when adding two dual numbers, we would like to append their name around the sign plus (+). Maybe string-append, let try:

scheme@(guile-user)> (define name1 "x")
scheme@(guile-user)> (define name2 "y")
scheme@(guile-user)> (string-append "(" name1 "+" name2 ")")
$9 = "(x+y)"

Neat! Well, thus for implementing the addition of dual numbers, we only need to access to all the fields and construct a new dual number,

(define (dual-add z1 z2)
  (let ((name1  (dual-name z1))
        (value1 (dual-value z1))
        (deriv1 (dual-deriv z1))
        (name2  (dual-name z2))
        (value2 (dual-value z2))
        (deriv2 (dual-deriv z2)))
    (dual
     (string-append "(" name1 "+" name2 ")")
     (+ value1 value2)
     (+ deriv1 deriv2))))

Well, something does not appear smooth. It could be nice to easy the unwrap way.

Pattern matching is a powerful feature by functional programming languages. It allows to match an object against several patterns and extract the elements that make it up. It requires the module (ice-9 match). Let take some examples to fix the ideas. The procedure example returns,

scheme@(guile-user)> (example (list 'first 2 'third))
$10 = third
scheme@(guile-user)> (example (list 'first 'second 3))
$11 = 9
scheme@(guile-user)> (example (list 1 2 3))
$12 = 5
scheme@(guile-user)> (example (list 1 2 3 5))
$13 = (1 2 3 5)
scheme@(guile-user)> (example "that")
$14 = "that"

How do you implement it? Using many conditionals? Nope, using pattern matching! The syntax reads,

(define (example object)
  (match object
    (('first 2 something)               ;1
     something)
    ((_ 'second n)                      ;2
     (* n n))
    ((1 b c)                            ;3
     (+ b c))
    (x                                  ;4
     x)))

where object is extracted following patterns:

1. if object matches the pattern ('first 2 something) where something is the name of a local variable, then return the value of that variable,

2. if object matches the pattern (_ 'second number) where _ (underscore) means anything and number is the name of a local variable, then return the multiplication with this number. It would be more correct to also check, when unwrapping that the third element is a number (predicate number?),

   ((_ 'second (? number? n))
    (* n n))
   

3. etc.

As you see, it is really powerful because on one hand, it extracts the elements and on the other hands it checks the pattern. Ok, another example. Lists are fundamental and primitive operations in Lisp are car and cdr which returns the first element of a list and the rest except the first.

scheme@(guile-user)> (car (list 'first 2 'third))
$15 = first
scheme@(guile-user)> (cdr (list 'first 2 'third))
$16 = (2 third)

The pattern matching allow ellipsis (...). It means that the previous pattern may be matched zero or more times in a list. Hence, a procedure using match returning the first element.

(define (car-match lst)
  (match lst
    ((head tail ...)
     head)))                            ;here

Instead of head (line ;here), try to replace by tail (and rename the procedure, e.g., cdr-match).

scheme@(guile-user)> (cdr-match (list 'first 2 'third))
$17 = (2 third)

So far, so good! Let implement the subtraction.

Once the pattern matching documentation careful read, instead of several access to the record field, now it appears straightforward to pattern match using the record. It reads,

(define (dual-sub z1 z2)
  (match z1
    (($ <dual> name1 value1 deriv1)
     (match z2
       (($ <dual> name2 value2 deriv2)
        (dual
         (string-append "(" name1 "-" name2 ")")
         (- value1 value2)
         (- deriv1 deriv2)))))))

Much more elegant than dual-add, isn’t it? The two remaining arithmetic operations (multiplication and division) are let as an exercise.

Since the exponential is an analytic function and we know its derivative (\((e^{u})^{\prime}=u^{\prime}e^{u}\)), the efficient implementation is to exploit such, for instance,

(define (dual-exp z)
  (match z
    (($ <dual> name value deriv)
     (let ((exp-value (exp value)))
       (dual (format #f "~~e^(~a)" (dual-name z))
             exp-value
             (* deriv exp-value))))))

Well, we would like to practise Guile, let implement the exponential using only the four arithmetic operations. The exponential function is defined by power series with the formula,

We will not compute until infinity, say only the first twenty terms. The first thinking is to iterate using a for-loop. However, it does not respect the functional style–instead recursion is the functional style. Well, we would say using plain words: if \(k\) is smaller than max-iter (say 20), then update and apply again (recursion-loop) this test, else return the result. Somehow, it would translate to,

(if (<= k max-iter)
    (loop (+1 k) expx xk kfac)
    (dual (format #f "~~e^(~a)" (dual-name x))
          (dual-value expx)
          (dual-deriv expx)))

where expx, xk and kfac are not well-defined yet. Knowing \(x^{k}\) or \(k!\), the update is provided by,

and thus, knowing previous summation (expx), we just compute the \(k\) -th term (term_k) and add it. Somehow, this loop procedure could be defined with,

(define (loop k expx xk kfac)
  (let* ((kk     (constant k))          ;convert integer to dual
         (xk     (dual-mul xk x))
         (kfac   (dual-mul kfac kk))
         (term_k (dual-div xk kfac))
         (expx   (dual-add expx term_k)))
    (if (<= k max-iter)
        (loop (1+ k) expx xk kfac)
        (dual (format #f "~~e^(~a)" (dual-name x))
              (dual-value expx)
              (dual-deriv expx)))))

And what remains is the initialization and/or initial call. It can be confusing at first but this loop procedure can be defined via let (see iteration mechanism).

(define (dual-exponential x)
  (let ((max-iter 20))
    (let loop ((k 1)
               (expx (constant 1))
               (xk   (constant 1))
               (kfac (constant 1)))
      (let* ((kk     (constant k))
             (xk     (dual-mul xk x))
             (kfac   (dual-mul kfac kk))
             (term_k (dual-div xk kfac))
             (expx   (dual-add expx term_k)))
        (if (<= k max-iter)
            (loop (1+ k) expx xk kfac)
            (dual (format #f "~~e^(~a)" (dual-name x))
                  (dual-value expx)
                  (dual-deriv expx)))))))

Confusing… yeah! Functional style appears unexpected at first; it is true for OCaml, Haskell or similar when coming from imperative style as Python or similar.

Because we have fun, let demo some other nice features.

In functional style, procedure (function) is a first-class citizen. Scheme makes it easy to compose them. Consider two procedures,

(define (f x)
  (+ x 2))

(define (g x)
  (* x 2))

and composing them reads,

(define (h x)
    (f (g x)))

or simply,

(define fog
  (compose f g))

And (fog 5) returns 12 as expected. Now, consider the new procedure f taking a value argument instead of plain 2.

(define (f x value)
  (+ x value))

And we would like to be able to vary this value in the composition fog. First, note that the composition does not match. Well, to understand correctly, we are going to write using good ol’ lambda,

(define fog
  (lambda (value)
    (lambda (x)
      (f (g x) value))))

Hum, what does it mean? Applying an argument (fog 2) returns a new procedure with one argument.

scheme@(guile-user)> (fog 2)
$18 = #<procedure 7fae6c192180 at <unknown port>:542:4 (x)>
scheme@(guile-user)> (define fog-2 (fog 2))
scheme@(guile-user)> (fog-2 5)
$19 = 12

The module (use-modules (srfi srfi-26)) provides a syntax for conveniently specializing selected parameters of a function; especially cut.

(define fog
  (lambda (value)
    (compose (cut f <> value) g)))

where the slot <> is filled by an argument from the created procedure.

Another example. From a number (potential floating-point), let convert to an integer. To do so, we convert the number as a string, then split the string against the dot character and last convert back the string to a number. The code reads,

(define number->integer
  (compose
     string->number car (cut string-split <> #\.) number->string))

And it works, (number->integer 1.2) returns 1.

As we said earlier, the list data structure plays a central role. Let manipulate it. It requires (use-modules (srfi srfi-1)).

The first essential component is map. Assuming a procedure and a list, the map procedure apply the procedure to each element of the list and returns a list. Roughly, naming f the procedure and considering the list (list 1 2 3), the result of (map f (list 1 2 3)) will be (list (f 1) (f 2) (f 3)).

The procedure map is built-in but now we have all the tools to have an implementation, for instance, renamed nap here,

(define (nap f lst)
  (match lst
    ((x xs ...)
     (cons (f x) (nap f xs)))
    (_ '())))

Give a look to this introduction about monads for some questions about composition.

Another essential component is fold (also termed reduce, accumulate, aggregate, compress, or inject). Assuming we would like to sum all the elements of the list (1 2 3 4 5), the structure of the operations reads,

  +
 / \
5   +
   / \
  4   +
     / \
    3   +
       / \
      2   +
         / \
        1   0

which is exactly what fold does. For instance, the code reads,

(define (sum-list lst)
  (fold (lambda (current previous)      ;reducer
          (+ current previous))
        0                               ;init
        lst))                           ;list

and the result is indeed 15.

scheme@(guile-user)> (sum-list (iota 5 1))
$20 = 15

where iota returns a sequence of 5 elements starting with 1.

It becomes straightforward to implement a procedure that returns (-3 -2 -1 0 1 2 3). It just reads,

(define (make-real-numbers n)
  (let ((positives (iota (1+ n))))
    (fold (lambda (element result)
            (cons (- element) result))
          (cdr positives)
          positives)))

where cons appends one element to an existing list. Based on that, we can also make dual numbers.

(define (make-variables n)
  (let* ((varname (lambda (i)
                   (string-append "x_" (number->string i)))))
    (map (lambda (i)
           (variable
            (varname i)
            i))
         (make-numbers n))))

and the outputs looks like,

scheme@(guile-user)> ,pp (make-variables 3)
$21 = (#<<dual> name: "x_-3" value: -3 deriv: 1.0>
 #<<dual> name: "x_-2" value: -2 deriv: 1.0>
 #<<dual> name: "x_-1" value: -1 deriv: 1.0>
 #<<dual> name: "x_0" value: 0 deriv: 1.0>
 #<<dual> name: "x_1" value: 1 deriv: 1.0>
 #<<dual> name: "x_2" value: 2 deriv: 1.0>
 #<<dual> name: "x_3" value: 3 deriv: 1.0>)

Awesome! One last?

scheme@(guile-user)> (for-each (compose dual-show dual-exponential2) (make-variables 2))
"~e^(x_-2)"
Value: 1.3533528323660934E-1
Deriv: 1.353352832366504E-1
"~e^(x_-1)"
Value: 3.6787944117144245E-1
Deriv: 3.6787944117144245E-1
"~e^(x_0)"
Value: 1.0E+0
Deriv: 1.0E+0
"~e^(x_1)"
Value: 2.7182818284590455E+0
Deriv: 2.7182818284590455E+0
"~e^(x_2)"
Value: 7.389056098930645E+0
Deriv: 7.389056098930604E+0

with,

(define (dual-exponential2 x)
  (define result
    (let ((max-iter 20))
      (fold (lambda (k result)
              (match result
                ((xk kfac expx)
                 (let* ((kk     (constant k))
                        (xk     (dual-mul xk x))
                        (kfac   (dual-mul kfac kk))
                        (term_k (dual-div xk kfac)))
                   (list
                    xk
                    kfac
                    (dual-add expx term_k))))))
            (list
             (constant 1.)
             (constant 1.)
             (constant 1.))
            (iota (1+ max-iter) 1))))

  (match result
    ((_ _ ($ <dual> _ value deriv))
     (dual (format #f "~~e^(~a)" (dual-name x))
           value
           deriv))))