7 ms·

Automatic Differentiation in 38 lines of Haskell

154 points by ttesmer 4y ago

mountainriver 4y ago
Being purely functional makes this quite easy, still a beauty to see
- sterlind 4y ago
  it's not just the functionalness. you couldn't do it this way in Lisp/Scheme (I think?) because of the lack of multiple dispatch. If you did (f 'x) for instance, you'd end up with things like (* 2 'x) which would blow up, since Lisp would try to compute the answer instead giving you '(* 2 x) back.
  - cryptonector 4y ago
    The Common Lisp Object System has multiple dispatch.
    - patrec 4y ago
      But math operations are not generic functions. So you'd need to create your own version of all math functions. And that won't compose with independently developed libraries (unlike the Haskell example elsewhere in the thread).
  - jhgb 4y ago
    I'm sure that systems like ScmUtils wouldn't have problems with this.
green_on_black 4y ago
Curious: I don't see the `^` op defined, or is it translated inti `exp` im guessing?
- ayjtyjtyjgfjc 4y ago
  It is one of three exponentiation operators in standard Haskell, no translation involved. https://stackoverflow.com/a/6400630/14768587 https://stackoverflow.com/a/6400630/14768587
- ttesmer 4y ago
  Since the type was made an instance of the Num typeclass, any function that can be used with Num's, can now be used on the type (Dual d). As per the Prelude[1], ^ is part of the Num typeclass. Same thing for * for Floating[2]. The hyperbolic tangent can also be used without being explicitly coded, as it can be derived using cosh and sinh! EDIT: As for the differentiation, it works for ^ since it is just multiplication (https://hackage.haskell.org/package/base-4.17.0.0/docs/src/GHC.Real.html#%5E https://hackage.haskell.org/package/base-4.17.0.0/docs/src/G...) for which the derivative was defined using the product rule. [1]: https://hackage.haskell.org/package/base-4.17.0.0/docs/Prelude.html#v:-94- https://hackage.haskell.org/package/base-4.17.0.0/docs/Prelu... [2]: https://hackage.haskell.org/package/base-4.17.0.0/docs/Prelude.html#v:-42--42- https://hackage.haskell.org/package/base-4.17.0.0/docs/Prelu...
- tikhonj 4y ago
  The ^ operator is defined in Haskell's standard library for raising values of any numeric type to non-negative, integral powers. Conceptually a ^ n just expands to a * a * ... * a, but the actual code is a bit more complex[1] for performance reasons. The neat thing with this approach is that ^ works for any numeric type, including user-defined types like Dual in this example. Since the Dual type can handle calculating derivatives for *, it gets derivatives for ^ for free. [1]: https://hackage.haskell.org/package/base-4.17.0.0/docs/src/GHC.Real.html#%5E https://hackage.haskell.org/package/base-4.17.0.0/docs/src/G...
  - green_on_black 4y ago
    Ah, `^` being naturals-only makes more sense in terms of how it could work with no special logic!
    - 4y ago
sterlind 4y ago
This is an interesting approach. Haskell is not a symbolic language, but you take advantage of the abstractness of type parameters in function definitions to thread your implementation of "D x" through, and pattern match on that. It's a neat design pattern. I bet it'd work in Julia too.
- xiphias2 4y ago
  Yes, but Julia has both forward and backward differention implemented (backwards it's harder).
  - ttesmer 4y ago
    As I wrote in the Markdown file, there's also a usable package for Haskell called `ad` on Hackage. It has both forward and backward autodiff and prevents expression swell, among other things. This gist is just for illustration purposes.
  - warinukraine 4y ago
- version_five 4y ago
  For an example in julia, see Mike Innes tutorial: https://github.com/MikeInnes/diff-zoo https://github.com/MikeInnes/diff-zoo I'm only a beginner in Julia and not and AD expert, but I went through the exercise of porting this to python and found it very enlightening
- bmitc 4y ago
  > Haskell is not a symbolic language I'm not sure I understand what this means. What is a symbolic language that excludes languages like Haskell, F#, OCaml, etc.?
  - klipt 4y ago
    Perhaps they mean not homoiconic like lisp or forth.
melony 4y ago
Where's the vector Jacobian matrix?
- omnicognate 4y ago
  You don't need to explicitly form a Jacobian matrix to do AD. The key insight is that the matrix is just a representation of a linear function. In forward mode AD you evaluate the linear function the Jacobian represents (i.e. the derivative) at the same time as evaluating the function you are differentiating, usually without ever explicitly building the matrix. The derivative is a local linear approximation of the original function, so it composes the same way. This allows you to compute the parts of it as you compute the parts of the original function, and pass on the intermediate results alongside the intermediate results of the original calculation.
mghwaz 4y ago
How do you come up with code looking this good? I've been working in Haskell for 6 months professionally: currently - I wouldn't be able to come up with something like this in a day (or more).
- chrsig 4y ago
  generally this isn't something that you come up with in a day. it's distilled down and refined over time, you just see it when it reaches maturity. can't speak to the OPs process though, maybe they shat it out on a whim :)
chrsig 4y ago
one of the references is a talk by simon peyton jones on AD[0] I'm neither a math expert nor a haskell expert, but I happen to enjoy both. It's been a while since I've watched a SPJ lecture, and I'd forgotten how much I want him to explain everything. [0] https://www.youtube.com/watch?v=EPGqzkEZWyw https://www.youtube.com/watch?v=EPGqzkEZWyw
sordina 4y ago
One of my favourite Haskell "one-liners" is combining the AD package with Number.Symbolic: {-# LANGUAGE ImportQualifiedPost #-} module Module_1663406024_9206 where import Numeric.AD qualified as Ad import Data.Number.Symbolic qualified as Sym -- >>> f x = x^2 + 3 * x -- >>> Ad.diff f 1 -- >>> Ad.diff f (Sym.var "a") -- 5 -- 3+a+a -- >>> Ad.diff sin pi -- >>> Ad.diff sin (Sym.var "a") -- -1.0 -- cos a The package authors did not need to coordinate to make this possible which is pretty wild. Saw it first here: https://twitter.com/GabriellaG439/status/647601518871359489 https://twitter.com/GabriellaG439/status/647601518871359489 and https://www.reddit.com/r/haskell/comments/3r75hq/comment/cwmqphe/?utm_source=share&utm_medium=web2x&context=3 https://www.reddit.com/r/haskell/comments/3r75hq/comment/cwm...
- cryptonector 4y ago
  > The package authors did not need to coordinate to make this possible which is pretty wild. It works because `f` is polymorphic. The type of its `x` argument is not constrained in `f`'s definition, so you can plug in any `x` of any type you want provided that `x`'s type implements the methods used in `f`'s definition. With the `Dual` scheme you get to use as `x` a "dual" of `y` (`f x`, for some `f`) and `y'`, and then you get an `f` applied to that `x` where the actual `f` is parameterized by the actual `x`'s type, and so the methods called by `f` are those that apply to `x`'s type. So instead of the traditional numeric addition and multiplication, you'd get the "dual" addition and multiplication, and then everything "chains" through and you end up with `diff f x` being the `y'` in the dual of `y` and `y'` (you don't care about the `y`, just the `y'` because you want the `diff` -- the differential or derivative). It's brilliant.
  - magicalhippo 4y ago
    I still don't get how sin ends up as cos, without any coordination.
    - krastanov 4y ago
      Presumably the Ad package has a list of known derivatives. The Sym package now "automatically" uses it, without ever having to have known of it. The "coordination" is that they both use the "symbol" sin to refer to the idea of sine function.
deleted 4y ago
[deleted]
tbensky 4y ago
Fun exercise in Prolog too: https://www.codebymath.com/index.php/welcome/lesson/prolog-deriv https://www.codebymath.com/index.php/welcome/lesson/prolog-d....
quickthrower2 4y ago
Is the Float' type there to not pollute the Float type with the defined typeclasses? Or is there something I am missing?
- bradrn 4y ago
  I have no clue. It looks totally redundant to me too.
  - ttesmer 4y ago
    You're absolutely right, it is totally redundant. You could reduce the entire thing even more in size by doing only `D Float Float`, removing the `TypeSynonymInstances` pragma and removing the `VectorSpace`, instead using normal Float operations for the derivatives in the Num, Floating and Fractional instances. Further, you can remove the `diff` function entirely if you want, since it's doing something very simple; calling `f` using a `D` instead of a normal Float, Int, etc. So you would simply get the function applied to x and it's derivative by doing `(\x -> x^2) (D 2 1)`, setting the derivative to 1 and x=2. However, I think the Float' and `diff` etc. is at least a little helpful in understanding it. I got it from SPJ's talk, which I linked to in the file. Also, it makes it easier (e.g. in the case of `diff`) to later add onto the Autodiff, for example by implementing reverse mode, Jacobians, etc.
MrUssek 4y ago
To be clear this is a forward mode auto diff implementation, not reverse mode, as might be inferred by the reference to the SPJ talk, correct?
- sriku 4y ago
  Correct. Reverse mode seems much harder to express in Haskell..i was trying to understand AD some time ago and used Haskell for that (before running into Conal's paper). This way of writing forward AD was easy and it was awesome to see type inference and laziness help me with the understanding. At that time, I tried to code up reverse mode AD and failed to do it with comparable simplicity. If I may ... 1. First attempt - https://sriku.org/blog/2019/03/08/automatic-differentiation/ https://sriku.org/blog/2019/03/08/automatic-differentiation/ 2. Dual numbers and Taylor numbers - http://sriku.org/blog/2019/03/12/automatic-differentiation-dual-numbers-taylor-numbers/ http://sriku.org/blog/2019/03/12/automatic-differentiation-d... 3. Higher ranked beings - http://sriku.org/blog/2019/03/13/automatic-differentiation-higher-ranked-beings/ http://sriku.org/blog/2019/03/13/automatic-differentiation-h...
  - sriku 4y ago
    The struggle I had with reverse mode felt comparable to my struggle with coding up union find (for which mutation is crucial for its efficiency).
  - iamwil 4y ago
    Which paper by Conal?
    - sriku 4y ago
      https://arxiv.org/abs/1804.00746 https://arxiv.org/abs/1804.00746
      - iamwil 4y ago
        Did this end up helping you with forward or reverse auto diff?
cryptonector 4y ago
log (D u u') = D (log u) (scale (log u) u') That doesn't look like the derivative of the natural logarithm!
- pjscott 4y ago
  (It was fixed shortly after your comment was posted.)
bmitc 4y ago
Forward-mode automatic differentiation is always fun to see because of the power but simplicity of the method. Any language with pattern matching makes it almost trivial to implement. Although, I'm rarely interested in <such and such> in <n> lines of <language>. The more interesting things are overall conciseness with regards to the problem, the expressiveness, and the clarity that the code produces.
naasking 4y ago
For those less familiar with Haskell, here's a C# implementation of dual numbers comparable to this: https://github.com/naasking/AutoDiffSharp/blob/master/AutoDiffSharp/Dual.cs https://github.com/naasking/AutoDiffSharp/blob/master/AutoDi...
mhh__ 4y ago
It's not as pretty but forward mode AD isn't that many lines even in C

154 points by ttesmer 4y ago

mountainriver 4y ago

Being purely functional makes this quite easy, still a beauty to see

sterlind 4y ago
it's not just the functionalness. you couldn't do it this way in Lisp/Scheme (I think?) because of the lack of multiple dispatch. If you did (f 'x) for instance, you'd end up with things like (* 2 'x) which would blow up, since Lisp would try to compute the answer instead giving you '(* 2 x) back.
- cryptonector 4y ago
  The Common Lisp Object System has multiple dispatch.
  - patrec 4y ago
    But math operations are not generic functions. So you'd need to create your own version of all math functions. And that won't compose with independently developed libraries (unlike the Haskell example elsewhere in the thread).
- jhgb 4y ago
  I'm sure that systems like ScmUtils wouldn't have problems with this.

green_on_black 4y ago

Curious: I don't see the `^` op defined, or is it translated inti `exp` im guessing?

ayjtyjtyjgfjc 4y ago
It is one of three exponentiation operators in standard Haskell, no translation involved. https://stackoverflow.com/a/6400630/14768587 https://stackoverflow.com/a/6400630/14768587
ttesmer 4y ago
Since the type was made an instance of the Num typeclass, any function that can be used with Num's, can now be used on the type (Dual d). As per the Prelude[1], ^ is part of the Num typeclass. Same thing for * for Floating[2]. The hyperbolic tangent can also be used without being explicitly coded, as it can be derived using cosh and sinh! EDIT: As for the differentiation, it works for ^ since it is just multiplication (https://hackage.haskell.org/package/base-4.17.0.0/docs/src/GHC.Real.html#%5E https://hackage.haskell.org/package/base-4.17.0.0/docs/src/G...) for which the derivative was defined using the product rule. [1]: https://hackage.haskell.org/package/base-4.17.0.0/docs/Prelude.html#v:-94- https://hackage.haskell.org/package/base-4.17.0.0/docs/Prelu... [2]: https://hackage.haskell.org/package/base-4.17.0.0/docs/Prelude.html#v:-42--42- https://hackage.haskell.org/package/base-4.17.0.0/docs/Prelu...
tikhonj 4y ago
The ^ operator is defined in Haskell's standard library for raising values of any numeric type to non-negative, integral powers. Conceptually a ^ n just expands to a * a * ... * a, but the actual code is a bit more complex[1] for performance reasons. The neat thing with this approach is that ^ works for any numeric type, including user-defined types like Dual in this example. Since the Dual type can handle calculating derivatives for *, it gets derivatives for ^ for free. [1]: https://hackage.haskell.org/package/base-4.17.0.0/docs/src/GHC.Real.html#%5E https://hackage.haskell.org/package/base-4.17.0.0/docs/src/G...
- green_on_black 4y ago
  Ah, `^` being naturals-only makes more sense in terms of how it could work with no special logic!
  - 4y ago

sterlind 4y ago

This is an interesting approach. Haskell is not a symbolic language, but you take advantage of the abstractness of type parameters in function definitions to thread your implementation of "D x" through, and pattern match on that. It's a neat design pattern. I bet it'd work in Julia too.

xiphias2 4y ago
Yes, but Julia has both forward and backward differention implemented (backwards it's harder).
- ttesmer 4y ago
  As I wrote in the Markdown file, there's also a usable package for Haskell called `ad` on Hackage. It has both forward and backward autodiff and prevents expression swell, among other things. This gist is just for illustration purposes.
- warinukraine 4y ago
version_five 4y ago
For an example in julia, see Mike Innes tutorial: https://github.com/MikeInnes/diff-zoo https://github.com/MikeInnes/diff-zoo I'm only a beginner in Julia and not and AD expert, but I went through the exercise of porting this to python and found it very enlightening
bmitc 4y ago
> Haskell is not a symbolic language I'm not sure I understand what this means. What is a symbolic language that excludes languages like Haskell, F#, OCaml, etc.?
- klipt 4y ago
  Perhaps they mean not homoiconic like lisp or forth.

melony 4y ago

Where's the vector Jacobian matrix?

omnicognate 4y ago
You don't need to explicitly form a Jacobian matrix to do AD. The key insight is that the matrix is just a representation of a linear function. In forward mode AD you evaluate the linear function the Jacobian represents (i.e. the derivative) at the same time as evaluating the function you are differentiating, usually without ever explicitly building the matrix. The derivative is a local linear approximation of the original function, so it composes the same way. This allows you to compute the parts of it as you compute the parts of the original function, and pass on the intermediate results alongside the intermediate results of the original calculation.

mghwaz 4y ago

How do you come up with code looking this good? I've been working in Haskell for 6 months professionally: currently - I wouldn't be able to come up with something like this in a day (or more).

chrsig 4y ago
generally this isn't something that you come up with in a day. it's distilled down and refined over time, you just see it when it reaches maturity. can't speak to the OPs process though, maybe they shat it out on a whim :)

chrsig 4y ago

one of the references is a talk by simon peyton jones on AD[0] I'm neither a math expert nor a haskell expert, but I happen to enjoy both. It's been a while since I've watched a SPJ lecture, and I'd forgotten how much I want him to explain everything. [0] https://www.youtube.com/watch?v=EPGqzkEZWyw https://www.youtube.com/watch?v=EPGqzkEZWyw

sordina 4y ago

One of my favourite Haskell "one-liners" is combining the AD package with Number.Symbolic: {-# LANGUAGE ImportQualifiedPost #-} module Module_1663406024_9206 where import Numeric.AD qualified as Ad import Data.Number.Symbolic qualified as Sym -- >>> f x = x^2 + 3 * x -- >>> Ad.diff f 1 -- >>> Ad.diff f (Sym.var "a") -- 5 -- 3+a+a -- >>> Ad.diff sin pi -- >>> Ad.diff sin (Sym.var "a") -- -1.0 -- cos a The package authors did not need to coordinate to make this possible which is pretty wild. Saw it first here: https://twitter.com/GabriellaG439/status/647601518871359489 https://twitter.com/GabriellaG439/status/647601518871359489 and https://www.reddit.com/r/haskell/comments/3r75hq/comment/cwmqphe/?utm_source=share&utm_medium=web2x&context=3 https://www.reddit.com/r/haskell/comments/3r75hq/comment/cwm...

cryptonector 4y ago
> The package authors did not need to coordinate to make this possible which is pretty wild. It works because `f` is polymorphic. The type of its `x` argument is not constrained in `f`'s definition, so you can plug in any `x` of any type you want provided that `x`'s type implements the methods used in `f`'s definition. With the `Dual` scheme you get to use as `x` a "dual" of `y` (`f x`, for some `f`) and `y'`, and then you get an `f` applied to that `x` where the actual `f` is parameterized by the actual `x`'s type, and so the methods called by `f` are those that apply to `x`'s type. So instead of the traditional numeric addition and multiplication, you'd get the "dual" addition and multiplication, and then everything "chains" through and you end up with `diff f x` being the `y'` in the dual of `y` and `y'` (you don't care about the `y`, just the `y'` because you want the `diff` -- the differential or derivative). It's brilliant.
- magicalhippo 4y ago
  I still don't get how sin ends up as cos, without any coordination.
  - krastanov 4y ago
    Presumably the Ad package has a list of known derivatives. The Sym package now "automatically" uses it, without ever having to have known of it. The "coordination" is that they both use the "symbol" sin to refer to the idea of sine function.

deleted 4y ago

[deleted]

tbensky 4y ago

Fun exercise in Prolog too: https://www.codebymath.com/index.php/welcome/lesson/prolog-deriv https://www.codebymath.com/index.php/welcome/lesson/prolog-d....

quickthrower2 4y ago

Is the Float' type there to not pollute the Float type with the defined typeclasses? Or is there something I am missing?

bradrn 4y ago
I have no clue. It looks totally redundant to me too.
- ttesmer 4y ago
  You're absolutely right, it is totally redundant. You could reduce the entire thing even more in size by doing only `D Float Float`, removing the `TypeSynonymInstances` pragma and removing the `VectorSpace`, instead using normal Float operations for the derivatives in the Num, Floating and Fractional instances. Further, you can remove the `diff` function entirely if you want, since it's doing something very simple; calling `f` using a `D` instead of a normal Float, Int, etc. So you would simply get the function applied to x and it's derivative by doing `(\x -> x^2) (D 2 1)`, setting the derivative to 1 and x=2. However, I think the Float' and `diff` etc. is at least a little helpful in understanding it. I got it from SPJ's talk, which I linked to in the file. Also, it makes it easier (e.g. in the case of `diff`) to later add onto the Autodiff, for example by implementing reverse mode, Jacobians, etc.

MrUssek 4y ago

To be clear this is a forward mode auto diff implementation, not reverse mode, as might be inferred by the reference to the SPJ talk, correct?

sriku 4y ago
Correct. Reverse mode seems much harder to express in Haskell..i was trying to understand AD some time ago and used Haskell for that (before running into Conal's paper). This way of writing forward AD was easy and it was awesome to see type inference and laziness help me with the understanding. At that time, I tried to code up reverse mode AD and failed to do it with comparable simplicity. If I may ... 1. First attempt - https://sriku.org/blog/2019/03/08/automatic-differentiation/ https://sriku.org/blog/2019/03/08/automatic-differentiation/ 2. Dual numbers and Taylor numbers - http://sriku.org/blog/2019/03/12/automatic-differentiation-dual-numbers-taylor-numbers/ http://sriku.org/blog/2019/03/12/automatic-differentiation-d... 3. Higher ranked beings - http://sriku.org/blog/2019/03/13/automatic-differentiation-higher-ranked-beings/ http://sriku.org/blog/2019/03/13/automatic-differentiation-h...
- sriku 4y ago
  The struggle I had with reverse mode felt comparable to my struggle with coding up union find (for which mutation is crucial for its efficiency).
- iamwil 4y ago
  Which paper by Conal?
  - sriku 4y ago
    https://arxiv.org/abs/1804.00746 https://arxiv.org/abs/1804.00746
    - iamwil 4y ago
      Did this end up helping you with forward or reverse auto diff?

cryptonector 4y ago

log (D u u') = D (log u) (scale (log u) u') That doesn't look like the derivative of the natural logarithm!

pjscott 4y ago
(It was fixed shortly after your comment was posted.)

bmitc 4y ago

Forward-mode automatic differentiation is always fun to see because of the power but simplicity of the method. Any language with pattern matching makes it almost trivial to implement. Although, I'm rarely interested in <such and such> in <n> lines of <language>. The more interesting things are overall conciseness with regards to the problem, the expressiveness, and the clarity that the code produces.

naasking 4y ago

For those less familiar with Haskell, here's a C# implementation of dual numbers comparable to this: https://github.com/naasking/AutoDiffSharp/blob/master/AutoDiffSharp/Dual.cs https://github.com/naasking/AutoDiffSharp/blob/master/AutoDi...

mhh__ 4y ago

It's not as pretty but forward mode AD isn't that many lines even in C