Derivative, dual numbers and Guile

introduction to some basics of Guile

Note: Back in 2019, I wrote «Calcul de Dérivées via les nombres duaux» using the C language. The aim was to illustrate how computing helps about concrete engineering problems. The same example is reused here to show Scheme programming. The context of this post is the monthly Café Guix. Join the fun!

In this post, we are going to illustrate some GNU Guile concepts. The programming language is based on Scheme and our aim is to be as functional and declarative as possible.

The jungle of the Lisp family programming language can be confusing, at first. Let quickly summarized. Lisp, initially named as « LISt Processor », is easily recognisable by its fully parenthesized prefix notation syntax. The first Lisp had been developed in 1958 and the modern variants inherit many components, making them all cousins. Scheme, introduced in the 1970s, is one variant and then various implementations of this programming language had been released; to name some Racket, MIT/GNU Scheme or Chez Scheme and… GNU Guile. Although language specifications are written down, all these implementations have their own specificity.

First thing first, the heart of any Lisp is S-expression, defined by an atom or an expression of the form (x y …) where x and y are S-expression. Some examples for the usual:

• (list 1 2 "foo") is the list containing the number 1, 2 and the string "foo",
• (define x (list 1 2 "foo")) binds the list to the symbol named x,
• (lambda (a) (+ 1 a)) creates a procedure taking one argument (named a) and then returning the addition by one,
• (define f (lambda (a) (+ 1 a))) binds the procedure to the symbol named f,
• the call of the procedure is simply: (f 42).

(define f
(lambda (a)
(+ 1 a)))


it is possible to use the equivalent sugar,

(define (f a)
(+ 1 a))


We have enough… well we will introduce later what we need. Let’s go!

The problem we try to solve (elementary maths)

This part is not necessary for the programming part. It just helps to have some context about what we are going to implement.

Let a mathematical function $$f$$ from the real numbers to the real numbers,

\begin{equation*} \begin{array}{lclcl} f &:& \mathbb{R} &\longrightarrow& \mathbb{R}\\ & & x &\mapsto& y=f(x) \end{array} \end{equation*}

and the derivative, noted $$f^\prime$$, is defined by,

$$\label{org6623050} f^\prime(x) = \underset{h\rightarrow 0}{\lim} \frac{f(x+h) - f(x)}{h}.$$

Picking one value, say $$a$$, it is straightforward to compute the function for this value and the new computed value is $$b=f(a)$$. When the function is analytic, symbolic computations are possible using pen and paper or some computer algebra system. From this symbolic formula of $$f^\prime$$, the value of the derivative can be then computed to the value $$a$$. However, most of the time, $$f$$ is not analytic and then how to compute $$f^\prime(a)$$?

For sure, it is possible to reuse the definition \eqref{org6623050} and approximate it by finite difference (fixing a small enough value for $$h$$). The main issue is the choice of such value $$h$$ and the accuracy strongly depends on this choice.

Differentiation rules

We learn in high school the differentiation rules which reads,

$$\label{orgb351b7d} \begin{array}{lcl} (f+g)^\prime &=& f^\prime + g^\prime\\ (f*g)^\prime &=& f^\prime * g + f* g^\prime\\ (f-g)^\prime &=& f^\prime - g^\prime\\ \displaystyle \left(\frac{f}{g}\right)^\prime &=& \displaystyle \frac{f^\prime * g - f* g^\prime}{g^2}. \end{array}$$

Nothing fancy!

Dual numbers

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,

\begin{equation*} x = (v, d) \end{equation*}

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

$$\label{orgaf50173} \begin{array}{lcll} \texttt{add}(x, y) &=& (v+u ,& d+p) \\ \texttt{mul}(x, y) &=& (u*v ,& d*v ~+~ u*p) \\ \texttt{sub}(x, y) &=& (u-v ,& d-p) \\ \texttt{div}(x, y) &=& (\displaystyle \frac{u}{v} ,& \displaystyle \frac{d*v ~-~ u*p}{p^2}). \end{array}$$

Wait! It looks really similar, not to 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.

Example

Let the function $$f(x)=x^{2}$$ and we are able to compute the derivative: $$f^{\prime}(x)=2x$$. Let the dual number $$z=(a, 1.)$$. Now, we substitute $$z$$ in $$f$$ applying the arithmetic operation of dual number,

\begin{equation*} \texttt{mul}(z,z) = (a*a, 1.*a + a*1.) = (a*a, 2*a). \end{equation*}

Tadam! Neat, isn’t it? Ok, we are to implement these dual numbers in Guile.

Dual number implemented with Guile

Record

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


Local variable and optional argument

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


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, the addition is thus only 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 have a way to unwrap. Pattern matching 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 number)                 ;2
(* number number))
((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.
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


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))
(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. Implementation of dual number exponential Since the exponential is an analytic function and we know its derivative ($$(e^{u})^{\prime}=u^{\prime}e^{u}$$), the correct 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,

\begin{equation*} e^{x} = \sum_{k=0}^{\infty} \frac{x^{k}}{k!} = 1 + x + \frac{x^{2}}{2} + \frac{x^{3}}{6} + \frac{x^{4}}{24} + \dots \end{equation*}

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,

\begin{equation*} \begin{array}{lcl|l} x^{k+1} &\leftarrow& x^{k} * x & x^{k}~\textrm{reads}~\texttt{xk}\\ (k+1)! &\leftarrow& k! * (k + 1) & k!~\textrm{reads}~\texttt{kfac} \end{array} \end{equation*}

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

Extra

High-order function

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


(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.

Lists: map, fold

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)))
(_ '())))


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


Take the pill

Join the fun, join Guix!

(last update: 2023-03-03 Fri 18:43)