Symbolic Differentiation in Pure
Johannes
in order to get familiar with Pure and to wrap my head around term
rewriting, I tried to implement symbolic differentiation in Pure as
an exercise. I am aware that there is a simple recipe in
http://pure-lang.googlecode.com/hg/pure/examples/symbolic.pure
, but the chain rule seems to be missing there. When implementing this
rule, the problem arises that we have to distinguish between variables
and constants. As an ad-hoc solution, I introduced a kind of
"Constructor" for constants which I call C, so e.g. the number 5 is
represented by (C 5). So I come up with the following rules:
using math ;
a * C 0 = C 0 ;
C 0 * a = C 0 ;
C 0 + a = a ;
a + C 0 = a ;
a * C 1 = a ;
C 1 * a = a ;
a^(C 1) = a ;
a^(C 0) = (C 1) ;
C a + C b = C (a + b) ;
C a * C b = C (a * b) ;
C a + c + C b = C (a + b) + c; /* rule 11 */
C a * c * C b = C (a * b) * c;
C a * (C b * c) = C (a * b) * c ;
c * C a * C b = c * C (a * b) ;
C a * c + C b * c = C (a + b) * c ;
c * C a + c * C b = C (a + b) * c ;
C a * ((C b * c) * d) = C (a * b) * c * d ;
C a * c + d + C b * c = C (a + b) * c + d ; /* rule 18 */
/* second argument of Deriv is the expression to be differentiated,
the first the differentiation variable */
Deriv x x = C 1 ;
Deriv x (a + b) = (Deriv x a) + (Deriv x b) ;
Deriv x (a * b) = b * (Deriv x a) + a * (Deriv x b) ;
Deriv x (sin y) = (cos y) * (Deriv x y) ;
Deriv x (cos y) = (-sin y) * (Deriv x y) ;
Deriv x (exp y) = (exp y) * (Deriv x y) ;
Deriv x (y ^ (C n)) =
if n == 0
then C 0
else (C n) * (y ^ (C (n-1))) * (Deriv x y) ;
Deriv x (C a) = C 0 ;
This works surprisingly well, for instance
> (Deriv x ((C 3)*x^(C 3)+(C 2)*(cos x)^(C 4)+(C 5)*x^(C 3))) ;
C 24*x^C 2+C 8*cos x^C 3*(-sin x)
or
> (Deriv x ((C 5) * (exp ((C 4) * cos (x^(C 2)))))) ;
C 5*(exp (C4*cos (x^C 2))*(C 4*((-sin (x^C 2))*(C 2*x))))
, so IMHO symbolic differentiation could be a nice example of the
capabilities of Pure and term rewriting! However:
1) The representation of constants looks ugly, of course. I tried to
replace the constructor function C by type tags :: double, but I
failed. Could this be done in any way?
2) My rules 11-18 are humbug, of course; they work in simplifying the
resulting expressions for the special examples I stated above, but
fail to simplify more complicated expressions. Is there a systematic
way for simplifying?
Best regards
Johannes Engels
Johannes
With
Deriv x (f y) = (Deriv y (f y)) * (Deriv x y) ;
Deriv x (sin x) = (cos x) ;
Deriv x (cos x) = (-sin x) ;
Deriv x (exp x) = (exp x) ;
Deriv x (x ^ (C n)) = if n == 0
then C 0
else (C n) * (x ^ (C (n-1))) ;
I get for my second example:
> Deriv x ((C 5)*(exp((C 4)*cos(x^(C 2))))) ;
C 5*Deriv x (exp (C 4*cos (x^C 2)))
instead of
that I have to implement the chain rule for each individual
differentiation rule as stated yesterday:
Deriv x (sin y) = (cos y) * (Deriv x y) ;
Deriv x (cos y) = (-sin y) * (Deriv x y) ;
Deriv x (exp y) = (exp y) * (Deriv x y) ;
Deriv x (y ^ (C n)) =
if n == 0
then C 0
else (C n) * (y ^ (C (n-1))) * (Deriv x y) ;
Cheers Johannes
John Cowan
> PS: Surprisingly, further abstraction of the chain rule does not work.
> With
>
> Deriv x (f y) = (Deriv y (f y)) * (Deriv x y) ;
"f" is in head position here, and is assumed to be literal.
--
Note that nobody these days would clamor for fundamental laws John Cowan
of *the theory of kangaroos*, showing why pseudo-kangaroos are co...@ccil.org
physically, logically, metaphysically impossible. http://www.ccil.org/~cowan
Kangaroos are wonderful, but not *that* wonderful. --Dan Dennett on zombies
Albert Graef
> 1) The representation of constants looks ugly, of course. I tried to
> replace the constructor function C by type tags :: double, but I
> failed. Could this be done in any way?
Maybe that failed because the constants were actually ints? :) There's a
generic number type which can be used to match anything that looks like
a number. (Pure has most of the usual numeric tower, see here for a list
of numeric types defined by the prelude:
http://docs.pure-lang.googlecode.com/hg/purelib.html#prelude-types.)
> 2) My rules 11-18 are humbug, of course; they work in simplifying the
> resulting expressions for the special examples I stated above, but
> fail to simplify more complicated expressions. Is there a systematic
> way for simplifying?
Simplifying arithmetic expressions can be surprisingly complex. It
depends on which kind of normal form you actually want to have. If terms
should be in sum of products form, you'll probably want to do some kind
of lexicographic sorting so that you end up with something like
c1*x11*...*x1n + c2*x21*...*x2m where the ci are constants and each
variable tuple (xi1,...,xin) only occurs once. This cannot be done
efficiently with simple toplevel rules for associativity, commutativity
and distributivity alone. You'll probably need to work with local rule
sets to get that kind of thing to work. See
http://docs.pure-lang.googlecode.com/hg/purelib.html#reduce.
Albert
--
Dr. Albert Gr"af
Dept. of Music-Informatics, University of Mainz, Germany
Email: Dr.G...@t-online.de, a...@muwiinfa.geschichte.uni-mainz.de
WWW: http://www.musikinformatik.uni-mainz.de/ag
Albert Graef
> Johannes scripsit:
>
>> PS: Surprisingly, further abstraction of the chain rule does not work.
>> With
>>
>> Deriv x (f y) = (Deriv y (f y)) * (Deriv x y) ;
>
> "f" is in head position here, and is assumed to be literal.
I should maybe point out that the right way to write this rule is:
Deriv x (f@_ y) = (Deriv y (f y)) * (Deriv x y) ;
See
http://docs.pure-lang.googlecode.com/hg/pure.html#head-function-pitfalls.
Johannes
On 29 Feb., 23:40, Albert Graef <Dr.Gr...@t-online.de> wrote:
> Maybe that failed because the constants were actually ints? :) There's a
> generic number type which can be used to match anything that looks like
> a number. (Pure has most of the usual numeric tower, see here for a list
> of numeric types defined by the prelude:http://docs.pure-lang.googlecode.com/hg/purelib.html#prelude-types.)
>
following code:
using math ;
Deriv x (a :: number) = 0 ;
Deriv x x = 1 ;
Deriv x (a * b) = b * (Deriv x a) + a * (Deriv x b) ;
Deriv x (a / b) = (b * (Deriv x a) - a * (Deriv x b)) / b^2 ;
/* quotient rule */
Deriv x (sin y) = (cos y) * (Deriv x y) ;
Deriv x (cos y) = (-sin y) * (Deriv x y) ;
Deriv x (exp y) = (exp y) * (Deriv x y) ;
if n == 0
then 0
else n * (y ^ (n-1)) * (Deriv x y) ;
/* some (crude) simplification rules: */
a * 0 = 0 ;
0 * a = 0 ;
a * 1 = a ;
a^0 = 1;
a :: number + c + b :: number = (a + b) + c;
a :: number * c * b :: number = (a * b) * c;
a :: number * (b :: number * c) = (a * b) * c ;
c * a :: number * b :: number = (a * b) * c ;
a :: number * c + b :: number * c = (a + b) * c ;
c * a :: number + c * b :: number = (a + b) * c ;
a :: number * ((b :: number * c) * d) = (a * b) * c * d ;
a :: number * c + d + b :: number * c = (a + b) * c + d ;
The implementation of the "generic" chain rule according to
as the chain rule can easily be implemented in the individual
differentiation rules as given above. The code above, which I find
amazingly simple, seems to cover already a lot of even nested
expressions, despite the crude simplification rules:
> Deriv x (3*(exp (4*(cos x)^2))) ;
3*(exp (4*cos x^2)*(8*cos x*(-sin x)))
>Deriv x (3*x^3 + 2*(cos x)^4 + 5*x^3) ;
24*x^2+8*cos x^3*(-sin x)
Further suggestions are welcome, of course.
Now I am really excited about Pure! Thank you very much for your help!
Best regards
Johannes Engels
Albert Graef
> The implementation of the "generic" chain rule according to
> Deriv x (f@_ y) = (Deriv y (f y)) * (Deriv x y) ;
> still failed (I got segmentation faults).
I haven't really reviewed your code in detail, but I suspect that the
above rule loops on 'Derive u (v w)', so you're getting a stack
overflow. You'll have to do more checking on the y argument to ensure
that the chain rule is really applicable.
Also, you might want to check
http://docs.pure-lang.googlecode.com/hg/pure.html#stack-size-and-tail-recursion
to see how you can set the PURE_STACK environment variable to get
orderly exceptions in case of a stack overflow. And then you might want
to read up on how to run your script with the symbolic debugger to see
where it goes wrong. ;-)
Grigory Sarnitskiy
>
> in order to get familiar with Pure and to wrap my head around term
> rewriting, I tried to implement symbolic differentiation in Pure as
> an exercise
thing. The exercise seems to be the natural point to start with :-)
But for now, I can only make a small remark. Once I've implemented
symbolic differentiation with Haskell using pattern matching with its
type constructors. However I didn't introduce 'x' and actually did the
whole thing in some sense pointfree. That is instead of 'x' i used id,
so
diff id = const 1
diff (const _) = const 0
diff (f + g) = diff f + diff g
diff (sin f) = (diff f) * (cos f)
I don't know yet whether it would be possible in pure, hope I'll be
able to find it out.
Grigory Sarnitskiy
diff id = cnst 1;
diff sin = cos;
to do what I want to? That is how to get rid of 'warning: rule never
reduced: diff sin = cos;'?
John Cowan
> diff id = cnst 1;
> diff sin = cos;
Pure distinguishes between functions/constructors (same thing) and
variables by whether they appear in the head of the expression or not.
So what you have written is the exact equivalent of:
diff x = cnst 1;
diff x = cos;
You can see why the second rule can never be reached.
The only way to manage a point-free style in Pure is to declare 'sin'
as a nonfix operator. But overall I think it is better to write "diff
sin x = cos x" and "diff x::number = 1".
When you try to write language A in the style of language B, you often
have problems like these. Just as with human languages, it is necessary
to learn to think in the second language to use it fluently.
--
John Cowan http://www.ccil.org/~cowan co...@ccil.org
"After all, would you consider a man without honor wealthy, even if his
Dinar laid end to end would reach from here to the Temple of Toplat?"
"No, I wouldn't", the beggar replied. "Why is that?" the Master asked.
"A Dinar doesn't go very far these days, Master. --Kehlog Albran
Besides, the Temple of Toplat is across the street." The Profit
Albert Graef
As John already pointed out, 'id' and 'sin' are interpreted as variables
there, as they are not the head of an application. You can fix that by
declaring them as "nonfix" symbols:
nonfix id sin;
But note that this changes the lexical meaning of those symbols so that
they are *always* interpreted as literals. To avoid that, you can either
use an explicit comparison:
diff f = cnst 1 if f === id;
= cos if f === sin;
Or you can "escape" the symbol by using a namespace qualifier:
diff ::id = cnst 1;
diff ::sin = cos;
See http://docs.pure-lang.googlecode.com/hg/pure.html#as-patterns (have
a look at the second part of that section which begins with "Sometimes
you may also run into the complementary problem, i.e., to match a
function argument against a given function.").
BTW, you may want to use cst instead of cnst. cst is defined in the
prelude as cst x y = x. See
http://docs.pure-lang.googlecode.com/hg/purelib.html#basic-combinators
Grigory Sarnitskiy
diff (f+g) = diff f + diff g;
diff (f*g) = f * diff g + g * diff f;
diff (f.g) = diff(f) * diff(g);
diff arg = cst 0 if arg === cst;
= cst 1 if arg === id;
= cos if arg === sin;
It works ok to a certain extent:
> diff $ id + id;
cst 1+cst 1
> diff $ id * sin + id;
id*cos+sin*cst 1+cst 1
> diff $ sin.sin;
cos*cos
But now I want to be able to apply the results as "functions":
> (diff $ sin) 1;
0.54030230586814
> (diff $ id) 12;
1
Obviously I can't right now do
> (diff $ id + sin) pi;
(cst 1+cos) 3.14159265358979
since I haven't defined + and * (pointwisely) for functions. And here
are two interconnected questions:
What is the difference between lambdas and equations? Is f x = sin x
the same as f = \x -> sin x ?
What is the best way to implement pointwise * and + for functions? so
I could have
>(cos + cos) 0;
2
? Shall I introduce types?
To John Cowan:
have problems like these. Just as with human languages, it is
necessary
to learn to think in the second language to use it fluently.
language I use to think about math. And I consider differentiation to
be a linear operator on functions, so points aren't fit well to this
point of view.
Ultimately I've decided to take pure because Haskell was a bit in
conflict with that 'language B'. Its notion of functions is to close
to the usual set-theoretic one. They are just mappings, the only thing
one can do with them is to apply to an argument. Contrary, I wanted a
bit of (extremely) old-fashioned 'symbolic' approach to them as
constructive definitions, which can be operated and rewritten.
Grigory Sarnitsky
using math;
f + g = f + g with (f+g) x = f x + g x end;
f * g = f * g with (f*g) x = f x * g x end;
f - g = f - g with (f-g) x = f x - g x end;
diff (f+g) = diff f + diff g;
diff (f*g) = f * diff g + g * diff f;
diff (f.g) = diff(f) * diff(g);
diff arg = cst 1 if arg === id;
= cos if arg === sin;
which has it merits:
> diff $ id + id;
cst 1+cst 1
cst 1+(id*cos+sin*cst 1)
however, now I can't fully apply + or *:
> (diff $ id + id) 1;
1+1
> (diff $ id + id*sin) 1;
1+(1*0.54030230586814+0.841470984807897*1)
By now I couldn't utilize + and * as
f + g = \x -> f x + g x;
f * g = \x -> f x * g x;
I still have a vague perception of the place of lambdas in term rewriting.
Albert Graef
> using math;
>
> f + g = f + g with (f+g) x = f x + g x end;
Hey, this trick is really neat, I wouldn't have thought of that! :) The
only problem is that the '+' on the rhs of the local definition is not
the same '+' as on the lhs. (Well, it's the same symbol, but it's bound
to the local '+' function there.) But that's easy to fix, you just need
to invoke the global '+' in the body of the local function:
f + g = f + g with (f+g) x = f x ::+ g x end;
NB: I'm not sure whether you need this here, but it might also be a good
idea to check that f and g are actually function objects which accept an
argument:
f + g = f + g with (f+g) x = f x ::+ g x end if nargs f > 0 && nargs g > 0;
If you do the same for the other operators, the rules will work just fine:
> (diff $ id + id) 1;
2
> (diff $ id + id*sin) 1;
2.38177329067604
> By now I couldn't utilize + and * as
>
> f + g = \x -> f x + g x;
> f * g = \x -> f x * g x;
No, this won't work because a lambda produces an anonymous function, and
there's no way you could match that in 'diff'.
> I still have a vague perception of the place of lambdas in term rewriting.
Well, lambdas don't exist in standard (first-order) term rewriting, they
are a Pure extension. And an evaluated lambda term can't be rewritten in
Pure, because it's just an atomic object (basically a function pointer)
as far as the rewriting machinery is concerned.
Pure only allows you to rewrite a lambda term while it's still
unevaluated, i.e., at macro expansion time. You might thus employ some
macro rules to wrap up the lambda in a way that would allow 'diff' to
dissect it, but that would make things much more complicated. I think
that your trick with the local '+' function is much nicer!
Albert Graef
> What is the difference between lambdas and equations? Is f x = sin x
> the same as f = \x -> sin x ?
Not quite. As you already found out, f evaluates to an anonymous
function in the latter case. Both will yield the same results on all
arguments, though.
Grigory Sarnitsky
f + g = f + g with (f+g) x = f x ::+ g x end if nargs f > 0 && nargs g > 0;
yields
> ((-id) + id) 1;
-1+1
That '-' denoting negate seems to be tricky in some languages. I wonder how things are in Pure.
The second hindrance is more substantial. Namely, my diff is ok for
> f1 = sin.cos + id;
> diff f1;
cos.cos*(-sin)+cst 1
but no wonder it is useless in the next case:
> f2 x = sin(cos x) + x;
> diff f2;
diff f2
That is in Pure I can pattern match expressions, but not rewriting rules. Maybe I'm stating the problem in wrong terms, but hope you can understand the right formulation from the examples. Is there a way to match rules?
Or maybe is there a way to have one-to-one correspondence between rules and expressions?
Maybe there were people who come from Haskell seeking for "pattern matching against lambdas" or for "(->) being an algebraic data type". If such, point please to the discussions.
Albert Graef
> f + g = f + g with (f+g) x = f x ::+ g x end if nargs f > 0 && nargs g > 0;
> -f = -f with (-f) x = neg (f x) end if nargs f > 0;
-x is a synonym for neg x, so you're actually defining a local neg
function there. Use ::neg on the rhs of the local definition and it will
work.
Albert Graef
> > diff f2;
> diff f2
You'll need reflection to deal with such cases. There are ways to
determine the rewriting rules defining a function:
> f2 x = sin(cos x) + x;
> get_fundef f2;
[f2 x-->sin (cos x)+x]
See http://docs.pure-lang.googlecode.com/hg/pure.html#reflection
At least in simple cases that would allow you to derive a point-free
form to which diff would be directly applicable. Of course in the
general case, if f2 is completely arbitrary, a closed-form differential
may not even exist.
Johannes
On 8 Apr., 17:14, Grigory Sarnitsky <yrog...@gmail.com> wrote:
>
> The second hindrance is more substantial. Namely, my diff is ok for
>
> > f1 = sin.cos + id;
> > diff f1;
>
> cos.cos*(-sin)+cst 1
>
> but no wonder it is useless in the next case:
>
> > f2 x = sin(cos x) + x;
> > diff f2;
>
> diff f2
>
> That is in Pure I can pattern match expressions, but not rewriting rules.
>
point-free style. With the code I posted on March, 2nd I get
> Deriv x (sin (cos x) + x);
cos (cos x)*(-sin x) +1
Best regards
Johannes
Michael Radford
definition and outputs parallel rules for the function's derivative.
But, the point-free style also breaks down with more than one
argument, whereas with explicit arguments, you can trivially
generalize to partial derivatives.
Mike
> --
> You received this message because you are subscribed to the Google Groups "pure-lang" group.
> To post to this group, send email to pure...@googlegroups.com.
> To unsubscribe from this group, send email to pure-lang+...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/pure-lang?hl=en.
>
Johannes Engels
rule (discussed on March, 2)
Deriv x (x ^ n) = if n == 0
then 0
else n * (x ^ (n-1)) ;
But e.g. ^ is actually a two-place function, i.e. for an expression
such as (2*x - 1)^3 the chain rule above fails, because "(f y)" is in
this case interpreted as a curried function. That means, I have to
modify the rule for ^.
As far as I see the following works:
using math ;
a * 0 = 0 ;
0 * a = 0 ;
a - 0 = a ;
0 - a = -a ;
0 ^ a = 0 ;
a ^ 1 = a ;
1 ^ a = 1 ;
Deriv x (sin x) = cos x ;
Deriv x (cos x) = -(sin x) ;
Deriv x (exp x) = exp x ;
/* ... add further elementary derivational rules here ... */
Deriv x (a :: number) = 0 ;
Deriv x (a + b) = (Deriv x a) + (Deriv x b) ; /* sum rule */
Deriv x (a - b) = (Deriv x a) - (Deriv x b) ;
Deriv x (a * b) = b * (Deriv x a) + a * (Deriv x b) ; /* product
rule*/
quotient rule */
(a^b * (ln a)) * (Deriv x b) ;
/* chain rules */
Deriv x (f@_ y z) = (Deriv y (f y z)) * (Deriv x y) +
(Deriv z (f y z)) * (Deriv x z) ;
Deriv x (f@_ y) = (Deriv y (f y)) * (Deriv x y) ;
> Deriv x (3*(x^3-1)^4 + 2*(cos (x^3))^4);
3*(4*(x^3-1)^3*(3*x^2))+2*(4*cos (x^3)^3*((-sin (x^3))*(3*x^2)))
The penultimate rule is probably redundant, but good to remember for
other two-place functions.
So, I have fun with Pure. Thanks again for all the useful hints!
Johannes