if {0} { [Arjen Markus] (3 may 2005) I came across the idea of ''automatic differentiation'' in a Fortran 90 program by [Simon Geard] that he was writing based on a Python example for the Language Shootout. While in Fortran 90 it is possible to define operators for so-called derived types, we can not do this in Tcl, at least not directly. The technique however is quite accessible and this page shows one way of doing it. First let me explain what the idea is: [Symbolic manipulation] is a kind of computer algebra that allows the computer to manipulate mathematical expressions, like rearranging the terms in the expression to make it more concise and easier to understand: a**2 + 2*a*b + b**2 ==> (a+b)**2 One kind of manipulation is to find the derivative of some function: f(x) = x**2 + 2*x df --(x) = 2*x + 2 dx The rules for that kind of manipulation are simple and can applied mechanically (in contrast: finding the primitive is a challenge that quite often leads to completely new classes of functions!) A very nice technique is to let the expression do the hard work: this is called ''automatic differentiation'' and here is one way to do it in Tcl: * First we define a bunch of procedures that allow us to express the function f above as a nesting of ordinary Tcl commands (these are easier to manipulate than the original expression): proc f {x} { return [+ [* \$x \$x] [* 2 \$x]] } * Then we use a list of two elements instead of plain numbers as arguments: set x [list 2.0 1.0] * The meaning of the first is: the ''value'' of x, the meaning of the second is: the value of the ''derivative'' of x with respect to the independent variable x. * Our procedures manipulate these two numbers, instead of the single value - this leads to the automatic evaluation of the function and its derivative. (Constants are represented as single numbers or as lists with 0.0 as the second element). So, the rules: f(x) = g(x) + h(x) ==> f'(x) = g'(x) + h'(x) f(x) = g(x) * h(x) ==> f'(x) = g'(x)*h(x) + g(x)*h'(x) etc. can be computed without us having to do more than specifying the way operators like + and * behave under differentiation.} namespace eval ::derivs { proc + {x y} { if { [llength \$x] == 1 } { set x [list \$x 0.0] } if { [llength \$y] == 1 } { set y [list \$y 0.0] } set v [expr {[lindex \$x 0]+[lindex \$y 0]}] set dv [expr {[lindex \$x 1]+[lindex \$y 1]}] return [list \$v \$dv] } proc * {x y} { if { [llength \$x] == 1 } { set x [list \$x 0.0] } if { [llength \$y] == 1 } { set y [list \$y 0.0] } set v [expr {[lindex \$x 0]*[lindex \$y 0]}] set dv [expr {[lindex \$x 1]*[lindex \$y 0] + [lindex \$x 0]*[lindex \$y 1]}] return [list \$v \$dv] } # - and / are left as an exercise ... # # Also define a few functions to make it more interesting # foreach {f df} {sin cos exp exp cos -sin log 1.0/ sqrt 0.5/sqrt} { proc \$f {x} [string map [list F \$f DF \$df] { if { [llength \$x] == 1 } { set x [list \$x 0.0] } set v [expr {F([lindex \$x 0])}] set dv [expr {[lindex \$x 1]*DF([lindex \$x 0])}] return [list \$v \$dv] }] } } # # A few tests: # Define two (global) functions # proc f {x} { return [* \$x [+ 2 \$x]] } proc g {x} { return [* 2 [* [sin \$x] [cos \$x]]] } namespace eval ::derivs { # # f'(x) = (2+x)+x = 2+2x # g'(x) = 2*cos(x)**2-2*sin(x)**2 = 2*cos(2x) # # "Import" the two functions # foreach p {f g} { proc \$p [info args ::\$p] [info body ::\$p] } # # Test the derivatives # puts "Parabola:" foreach x {0.0 0.5 1.0 1.5 2.0} { set xl [list \$x 1.0] ;# Very important!! puts "[f \$xl] -- [expr {2.0+2.0*\$x}]" } puts "sin(2x):" foreach x {0.0 0.5 1.0 1.5 2.0} { set xl [list \$x 1.0] puts "[g \$xl] -- [expr {2.0*cos(2.0*\$x)}]" } } if {0} { This can be extended to the evaluation of second derivatives if you like, or you can construct the expression that evaluates the derivative. Another, somewhat simpler, variation is to determine the ''size'' of the expression tree - something that is important with [genetic programming] techniques. Merely evaluate the above functions in a different namespace ... a neat trick [Donal Fellows] mentioned on comp.lang.tcl:} namespace eval ::treesize { foreach op {+ * - /} { # # Only the "arity" counts ... # proc \$op {x y} { set count 1 if { [llength \$x] == 1 } { incr count } else { incr count [lindex \$x 1] } if { [llength \$y] == 1 } { incr count } else { incr count [lindex \$y 1] } return [list COUNT \$count] } } # # Also define a few functions to make it more interesting # foreach f {sin exp cos} { proc \$f {x} { set count 1 if { [llength \$x] == 1 } { incr count } else { incr count [lindex \$x 1] } return [list COUNT \$count] } } } # # Reuse the functions above ... # namespace eval ::treesize { # # "Import" the two functions # foreach p {f g} { proc \$p [info args ::\$p] [info body ::\$p] } puts "Size of 'f': [f 1]" puts "Size of 'g': [g 1]" } if {0} { So, with the proper operators and mathematical functions we need only to define the actual expression once to automatically derive all kinds of mathematical and symbolic properties! Suppose you have an expression like: f(x) = ln( x + sqrt(1+x**2) ) and you need to compute the derivative, that is trivial with the code above. Manually verifying that its first derivative is: 1 f'(x) = ------------ sqrt(1+x**2) and then computing the value would take some time :)}