Recursion scheme for symbolic differentiation

Question

Following terminology from this excellent series, let's represent an expression such as (1 + x^2 - 3x)^3 by a Term Expr, where the data types are the following:

data Expr a =
    Var
  | Const Int
  | Plus a a
  | Mul a a
  | Pow a Int
  deriving (Functor, Show, Eq)

data Term f = In { out :: f (Term f) }

Is there a recursion scheme suitable for performing symbolic differentiation? I feel like it's almost a Futumorphism specialized to Term Expr, i.e. futu deriveFutu for an appropriate function deriveFutu:

data CoAttr f a  
  = Automatic a
  | Manual (f (CoAttr f a))

futu :: Functor f => (a -> f (CoAttr f a)) -> a -> Term f  
futu f = In <<< fmap worker <<< f where  
    worker (Automatic a) = futu f a
    worker (Manual g) = In (fmap worker g)

This looks pretty good, except that the underscored variables are Terms instead of CoAttrs:

deriveFutu :: Term Expr -> Expr (CoAttr Expr (Term Expr))
deriveFutu (In (Var))      = (Const 1)
deriveFutu (In (Const _))  = (Const 0)
deriveFutu (In (Plus x y)) = (Plus (Automatic x) (Automatic y))
deriveFutu (In (Mul x y))  = (Plus (Manual (Mul (Automatic x) (Manual _y)))
                                   (Manual (Mul (Manual _x) (Automatic y)))
                             )
deriveFutu (In (Pow x c))  = (Mul (Manual (Const c)) (Manual (Mul (Manual (Pow _x (c-1))) (Automatic x))))

The version without recursion schemes looks like this:

derive :: Term Expr -> Term Expr
derive (In (Var))      = In (Const 1)
derive (In (Const _))  = In (Const 0)
derive (In (Plus x y)) = In (Plus (derive x) (derive y))
derive (In (Mul x y))  = In (Plus (In (Mul (derive x) y)) (In (Mul x (derive y))))
derive (In (Pow x c))  = In (Mul (In (Const c)) (In (Mul (In (Pow x (c-1))) (derive x))))

As an extension to this question, is there a recursion scheme for differentiating and eliminating "empty" Exprs such as Plus (Const 0) x that arise as a result of differentiation -- in one pass over the data?

So your expressions always have one variable? – Willem Van Onsem May 22 '18 at 17:11 — Willem Van Onsem, May 22 '18 at 17:11
At most one, yes, if I understand your question correctly. – nnnmmm May 22 '18 at 17:13 — nnnmmm, May 22 '18 at 17:13

Roman Cheplyaka · Accepted Answer · 2018-05-22T18:32:57.510

Look at the differentiation rule for product:

(u v)' = u' v + v' u

What do you need to know to differentiate a product? You need to know the derivatives of the subterms (u', v'), as well as their values (u, v).

This is exactly what a paramorphism gives you.

para
  :: Functor f
  => (f (b, Term f) -> b)
  -> Term f -> b
para g (In a) = g $ (para g &&& id) <$> a

derivePara :: Term Expr -> Term Expr
derivePara = para $ In . \case
  Var -> Const 1
  Const _ -> Const 0
  Plus x y -> Plus (fst x) (fst y)
  Mul x y -> Plus
    (In $ Mul (fst x) (snd y))
    (In $ Mul (snd x) (fst y))
  Pow x c -> Mul
    (In (Const c))
    (In (Mul
      (In (Pow (snd x) (c-1)))
      (fst x)))

Inside the paramorphism, fst gives you access to the derivative of a subterm, while snd gives you the term itself.

As an extension to this question, is there a recursion scheme for differentiating and eliminating "empty" Exprs such as Plus (Const 0) x that arise as a result of differentiation -- in one pass over the data?

Yes, it's still a paramorphism. The easiest way to see this is to have smart constructors such as

plus :: Term Expr -> Term Expr -> Expr (Term Expr)
plus (In (Const 0)) (In x) = x
plus (In x) (In (Const 0)) = x
plus x y = Plus x y

and use them when defining the algebra. You could probably express this as some kind of para-cata fusion, too.

Yep. I came to the same conclusion, but then wasted time trying to find a way (within `recursion-schemes`) to get rid of the `In`s on the RHS. If there's a way, I couldn't find it. — dfeuer, May 22 '18 at 18:25
I think you reverse the pair compared to `recursion-schemes`, which might not be the best thing. I also think it might be better to pattern match on the pairs explicitly, with `x` and `dx` in whatever order. — dfeuer, May 22 '18 at 18:26
I think this still has a bug, it returns 3*(1+x²-3x)²*(1+x²-3x) for input (1+x²-3x)³ -- edit: Nevermind, I was using the `para` from the linked series which has its arguments reversed. — nnnmmm, May 22 '18 at 18:27

Recursion scheme for symbolic differentiation

1 Answers1