Safe forward-mode AD in Haskell?

Dear all, = Introduction = In [1] Siskind and Pearlmutter expose the danger of "perturbation confusion" in forward-mode automatic differentiation. In particular they state: "We discuss a potential problem with forward-mode AD common to many AD systems, including all attempts to integrate a forward-mode AD operator into Haskell." This literate Haskell message shows how, using a type system extension of the Glasgow Haskell Compiler[2] (GHC), we can statically guarantee that perturbation confusion does not occur. = Code = The below code closely follows the Haskell code in the appendices of [1]. It relies on GHC's arbitrary-rank polymorphism(?) extension to Haskell 98.
{-# OPTIONS_GHC -fglasgow-exts #-}
In our definition of the 'Bundle' data type we add the phantom type 's', which will be the key to disambiguating between different application of the 'd' operator.
data Bundle s a = Bundle a a
instance Num a => Show (Bundle s a) where showsPrec p (Bundle x x') = showsPrec p [x,x']
instance Num a => Eq (Bundle s a) where (Bundle x x') == (Bundle y y') = (x == y)
lift z = Bundle z 0
instance Num a => Num (Bundle s a) where (Bundle x x') + (Bundle y y') = Bundle (x + y) (x' + y') (Bundle x x') * (Bundle y y') = Bundle (x * y) (x * y' + x' * y) fromInteger z = lift (fromInteger z)
instance Fractional a => Fractional (Bundle s a) where fromRational z = lift (fromRational z)
We provide a type signature for 'd' where we existentially quantify the phantom type 's' to prevent mixing of bundles from different 'd' operators.
d :: Num a => (forall s. Bundle s a -> Bundle s a) -> a -> a d f x = let (Bundle y y') = f (Bundle x 1) in y'
The extential quantification makes the definition ] constant_one' x = d (\y -> x + y) 1 impossible since 'x' originates externally to the 'd' operator. GHC rejects the definition with a "Inferred type is less polymorphic than expected". In order to pass the compiler the definition must be changed to the corrected[1] version.
constant_one x = d (\y -> (lift x) + y) 1 should_be_one_a = d (\x -> x * (constant_one x)) 1 should_be_one_b = d (\x -> x * 1 ) 1
violation_of_referential_transparency = should_be_one_a /= should_be_one_b
= References = [1] http://www.bcl.hamilton.ie/~qobi/nesting/papers/ifl2005.pdf [2] http://haskell.org/ghc/

Oops, that went to soon.
My intention was to ask for comments from those who know better than
I. These corners of the type system are pretty unknown to me but I am
trying to learn from e.g. Oleg's array branding[1].
Anyway, please rip my exposition apart and expose my sure-to-be flawed
logic and terminology. ;)
Cheers,
Bjorn
[1] http://okmij.org/ftp/Haskell/eliminating-array-bound-check.lhs
On 5/8/07, Björn Buckwalter
Dear all,
= Introduction =
In [1] Siskind and Pearlmutter expose the danger of "perturbation confusion" in forward-mode automatic differentiation. In particular they state:
"We discuss a potential problem with forward-mode AD common to many AD systems, including all attempts to integrate a forward-mode AD operator into Haskell."
This literate Haskell message shows how, using a type system extension of the Glasgow Haskell Compiler[2] (GHC), we can statically guarantee that perturbation confusion does not occur.
= Code =
The below code closely follows the Haskell code in the appendices of [1]. It relies on GHC's arbitrary-rank polymorphism(?) extension to Haskell 98.
{-# OPTIONS_GHC -fglasgow-exts #-}
In our definition of the 'Bundle' data type we add the phantom type 's', which will be the key to disambiguating between different application of the 'd' operator.
data Bundle s a = Bundle a a
instance Num a => Show (Bundle s a) where showsPrec p (Bundle x x') = showsPrec p [x,x']
instance Num a => Eq (Bundle s a) where (Bundle x x') == (Bundle y y') = (x == y)
lift z = Bundle z 0
instance Num a => Num (Bundle s a) where (Bundle x x') + (Bundle y y') = Bundle (x + y) (x' + y') (Bundle x x') * (Bundle y y') = Bundle (x * y) (x * y' + x' * y) fromInteger z = lift (fromInteger z)
instance Fractional a => Fractional (Bundle s a) where fromRational z = lift (fromRational z)
We provide a type signature for 'd' where we existentially quantify the phantom type 's' to prevent mixing of bundles from different 'd' operators.
d :: Num a => (forall s. Bundle s a -> Bundle s a) -> a -> a d f x = let (Bundle y y') = f (Bundle x 1) in y'
The extential quantification makes the definition
] constant_one' x = d (\y -> x + y) 1
impossible since 'x' originates externally to the 'd' operator. GHC rejects the definition with a "Inferred type is less polymorphic than expected". In order to pass the compiler the definition must be changed to the corrected[1] version.
constant_one x = d (\y -> (lift x) + y) 1 should_be_one_a = d (\x -> x * (constant_one x)) 1 should_be_one_b = d (\x -> x * 1 ) 1
violation_of_referential_transparency = should_be_one_a /= should_be_one_b
= References =
[1] http://www.bcl.hamilton.ie/~qobi/nesting/papers/ifl2005.pdf [2] http://haskell.org/ghc/

On Tue, May 08, 2007 at 06:06:56PM -0400, Björn Buckwalter wrote:
We provide a type signature for 'd' where we existentially quantify the phantom type 's' to prevent mixing of bundles from different 'd' operators.
d :: Num a => (forall s. Bundle s a -> Bundle s a) -> a -> a d f x = let (Bundle y y') = f (Bundle x 1) in y'
Couldn't you do the same without branding if you simply made your function polymorphic: data Bundle a = Bundle a a d :: Num a => (forall b. Num b => b -> b) -> a -> a d f x = let (Bundle y y') = f (Bundle x 1) in y' and defined lift x = Bundle x 0 ? You obviously need an existential type, but I don't see that you need any branding. -- David Roundy Department of Physics Oregon State University

On Tue, May 08, 2007 at 05:59:48PM -0700, David Roundy wrote:
Couldn't you do the same without branding if you simply made your function polymorphic:
data Bundle a = Bundle a a
d :: Num a => (forall b. Num b => b -> b) -> a -> a d f x = let (Bundle y y') = f (Bundle x 1) in y'
and defined
lift x = Bundle x 0
I take it back. This lift won't work. You need a lift with the type lift :: (Num a, Num b) -> a -> b which would obviously need to be added to a class that would take the place of Num here... -- David Roundy Department of Physics Oregon State University

Okay, here's an alternate, unbranded approach:
data Bundle a = Bundle a a
instance Num a => Show (Bundle a) where showsPrec p (Bundle x x') = showsPrec p [x,x']
instance Num a => Eq (Bundle a) where (Bundle x x') == (Bundle y y') = (x == y)
instance Num a => Num (Bundle a) where (Bundle x x') + (Bundle y y') = Bundle (x + y) (x' + y') (Bundle x x') * (Bundle y y') = Bundle (x * y) (x * y' + x' * y) fromInteger z = Bundle (fromInteger z) 0
lift z = Bundle z 0
d :: Num a => (forall b. (Num b) => (a -> b) -> b -> b) -> a -> a d f x = let (Bundle y y') = f lift (Bundle x 1) in y'
The key change is in the type of d, which now accepts a polymorphic function on numbers, but passes in a "lift" function, which allows us to pass in higher-level variables. In one sense this function is ugly. In another sense, it's prettier, as you can now hide *all* of the Bundle data in a "differentiation" module.
constant_one x = d (\l y -> l x + y) 1
should_be_one_a = d (\_ x -> x * (constant_one x)) 1 should_be_one_b = d (\_ x -> x * 1 ) 1
violation_of_referential_transparency = should_be_one_a /= should_be_one_b
-- David Roundy Department of Physics Oregon State University
participants (2)
-
Björn Buckwalter
-
David Roundy