Tuesday, January 15, 2008

Deriving Miss Daisy

As I've mentioned previously, I'm slowly reading through SICP - it's an amazing book and reading it has been a humbling experience.

'Geeking' Out

I completely geeked out reading one example on functional derivatives.  (Take a math expression as a parameter and return the derivative as an expression)

It's in the first chapter, and they use this to implement Newton's Method -> to then implement a Sqrt function!

Approximate Derivative

After I picked my brain off the floor, I decided to have a crack at implementing derivatives in C# using lambda expressions:

derivative

Using the formula of a derivative (above), we can simply code an approximate derivative (below):

Func<double, double> Approximate(Func<double, double> f)
{
var h = 0.0001;
return a => ((f(a + h) - f(a)) / h);
}

Symbolic Derivation

Easy huh!  So for fun I decided to use linq expression trees to derive algebraically.

It's easy to get an expression tree in C# - simply use the Expression<> type.
Func<double, double> Derive(Expression<Func<double, double>> e)
{
var lambda = new Derivative().Eval(e) as LambdaExpression;
return (Func<double, double>) lambda.Compile();
}

I parse the expression tree using Matt Warren's ExpressionVisitor class - and using the basic rules of Derivation (Constant, Product and Quotient rules) we get:
class Derivative : ExpressionVisitor
{
public Expression Eval(Expression exp)
{
return this.Visit(Evaluator.PartialEval(exp));
}

protected override Expression VisitBinary(BinaryExpression b)
{
// Product Rule: (fg)' = f'g + fg'
if (b.NodeType == ExpressionType.Multiply)
{
return Expression.Add(
Expression.Multiply(Eval(b.Left), b.Right),
Expression.Multiply(b.Left, Eval(b.Right)));
}

// Quotient Rule: (f/g)' = (f'g - fg') / (g*g)
if (b.NodeType == ExpressionType.Divide)
{
return Expression.Divide(
Expression.Subtract(
Expression.Multiply(Eval(b.Left), b.Right),
Expression.Multiply(b.Left, Eval(b.Right))),
Expression.Multiply(b.Right, b.Right));
}

return base.VisitBinary(b);
}

// Parameter Derivation: f(x)' = 1
protected override Expression VisitParameter(ParameterExpression p)
{
return Expression.Constant(1.0);
}

// Constant Rule f(a)' = 0
protected override Expression VisitConstant(ConstantExpression c)
{
return Expression.Constant(0.0);
}
}

Wrapping it up

To quickly finish this long blog post, here's the usage code:
// f(x * x * x)' = 3 * x * x
var approx = Approximate(x => x * x * x);
var derive = Derive(x => x * x * x);

var y1 = approx(1); // 3.0000300001109532
var y2 = derive(1); // 3.0

And there it is!  Useful? maybe not.  Fun? Definitely!

4 comments:

nico said...

Hey nice example, and quite educational!
I assume that the Eval will derive recursively, is that so? When reading it at first glance it looks a bit like it would evaluate the function..

return Expression.Add(
Expression.Multiply(Eval(b.Left), b.Right),
Expression.Multiply(b.Left, Eval(b.Right)));

Wouldn't it be nice to have even some functional optimizer??
i.e. code that reduces the number of calculations in the tree, like
it matches a^2 + 2ab + b^2 (having 6 operations) and replaces it by (a+b)^2
i'll revisit your blog for any thoughts on it ;)

Luke Marshall said...

Thanks!

Yes, Eval will derive recursively.

Optimizer would be very cool. I actually started on one, but stopped after a few hours - it was getting quite complex!

nico said...

yes thats true, perhaps it helps to have a look on all those lisp coders' stuff around, they use "expression trees" all over :)
the fun things is (i've looked on one once) that of course thats all recursive and salted with lambda expressions, exactly what's all possible with the new extensions theres huge codebases on the net...
have a nice and refreshing weekend, anyhow

Tim said...

This reminds me of something I had to write in LISP in uni (and it melted my brain).

LOL. Nico's suggestion sounds like taking optimisation to the extreme... but bring it on! I'm sure a parallel expression tree could be constructed on the first visit and substituted for further iterations. I think the writer would have to wrap their head in a super-cooled cloth first, and have air-vents installed in the side of their head.