Re: [Haskell-cafe] Hit a wall with the type system

2007-11-30 Thread David Roundy
On Thu, Nov 29, 2007 at 04:25:43PM -0800, Dan Weston wrote:
 I must be missing something, because to me the contract seems to be much 
 simpler to express (than the Functor + Isomorphism route you seem to me 
 to be heading towards):
...
 diff f a = if dx == dx' then error Zero denom else dydx
   where a'   = addEpsilon  a
 dx   = a'   `takeAway`   a
 dx'  = a`takeAway`   a'
 dy   = f a' `takeAway` f a
 dydx = dy   `divide`   dx

But this is just a finite-difference derivative, which is generally
accurate to... well, a lot less than an analytic derivative.  e.g try
computing

diff cos 1e-30

you'll get a garbage answer, while AD or symbolic differentiation will give
you the correct answer, 1e-30.

So for derivatives you actually intend to use, it's much better to use AD
or symbolic differentiation, or just compute the derivatives by hand.
Anything but finite difference (except at a carefully examined check for
better derivatives).
-- 
David Roundy
Department of Physics
Oregon State University
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe


Re: [Haskell-cafe] Hit a wall with the type system

2007-11-29 Thread Henning Thielemann

On Wed, 28 Nov 2007, Chris Smith wrote:

  data AD a = AD a a deriving Eq
 
  instance Show a = Show (AD a) where
  show (AD x e) = show x ++  +  ++ show e ++  eps
 
  instance Num a = Num (AD a) where
  (AD x e) + (AD y f)   = AD (x + y) (e + f)
  (AD x e) - (AD y f)   = AD (x - y) (e - f)
  (AD x e) * (AD y f)   = AD (x * y) (e * y + x * f)
  negate (AD x e)   = AD (negate x)  (negate e)
  abs (AD 0 _)  = error not differentiable: |0|
  abs (AD x e)  = AD (abs x) (e * signum x)
  signum (AD 0 e)   = error not differentiable: signum(0)
  signum (AD x e)   = AD (signum x)  0
  fromInteger i = AD (fromInteger i) 0
 
  instance Fractional a = Fractional (AD a) where
  (AD x e) / (AD y f)   = AD (x / y) ((e * y - x * f) / (y * y))
  recip (AD x e)= AD (1 / x) ((-e) / (x * x))
  fromRational x= AD (fromRational x) 0
 
  instance Floating a = Floating (AD a) where
  pi= AD pi0
  exp (AD x e)  = AD (exp x)   (e * exp x)
  sqrt (AD x e) = AD (sqrt x)  (e / (2 * sqrt x))
  log (AD x e)  = AD (log x)   (e / x)
  (AD x e) ** (AD y f)  = AD (x ** y)  (e * y * (x ** (y-1)) +
f * (x ** y) * log x)
  sin (AD x e)  = AD (sin x)   (e * cos x)
  cos (AD x e)  = AD (cos x)   (-e * sin x)
  asin (AD x e) = AD (asin x)  (e / sqrt (1 - x ** 2))
  acos (AD x e) = AD (acos x)  (-e / sqrt (1 - x ** 2))
  atan (AD x e) = AD (atan x)  (e / (1 + x ** 2))
  sinh (AD x e) = AD (sinh x)  (e * cosh x)
  cosh (AD x e) = AD (cosh x)  (e * sinh x)
  asinh (AD x e)= AD (asinh x) (e / sqrt (x^2 + 1))
  acosh (AD x e)= AD (acosh x) (e / sqrt (x^2 - 1))
  atanh (AD x e)= AD (atanh x) (e / (1 - x^2))
 
  diffNum:: Num b= (forall a. Num a= a - a) - b 
  - b
  diffFractional :: Fractional b = (forall a. Fractional a = a - a) - b 
  - b
  diffFloating   :: Floating b   = (forall a. Floating a   = a - a) - b 
  - b
 
  diffNum f x= let AD y dy = f (AD x 1) in dy
  diffFractional f x = let AD y dy = f (AD x 1) in dy
  diffFloating f x   = let AD y dy = f (AD x 1) in dy

Why do the functions have different number types after differentiation?
I thought, that just 'diff'

diff :: (AD a - AD a) - (a - a)
diff f x   = let AD y dy = f (AD x 1) in dy

must work. What you do, looks like numbers with errors, but I suspect you
are right that 'automatic differentiation' is the term used for those
kinds of computations.
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe


Re: [Haskell-cafe] Hit a wall with the type system

2007-11-29 Thread Henning Thielemann

On Thu, 29 Nov 2007, Henning Thielemann wrote:

 On Wed, 28 Nov 2007, Chris Smith wrote:

   diffNum:: Num b= (forall a. Num a= a - a) - b 
   - b
   diffFractional :: Fractional b = (forall a. Fractional a = a - a) - b 
   - b
   diffFloating   :: Floating b   = (forall a. Floating a   = a - a) - b 
   - b
  
   diffNum f x= let AD y dy = f (AD x 1) in dy
   diffFractional f x = let AD y dy = f (AD x 1) in dy
   diffFloating f x   = let AD y dy = f (AD x 1) in dy

 Why do the functions have different number types after differentiation?
 I thought, that just 'diff'

 diff :: (AD a - AD a) - (a - a)
 diff f x   = let AD y dy = f (AD x 1) in dy

 must work. What you do, looks like numbers with errors, but I suspect you
 are right that 'automatic differentiation' is the term used for those
 kinds of computations.

I like to add that I have written some code for working with Taylor
expansions of functions, that is, with all derivatives of a function. You
can get the derivatives of composed functions by composing the Taylor
series of the elementary functions:
  http://darcs.haskell.org/htam/src/PowerSeries/Taylor.hs

E.g.

*PowerSeries.Taylor take 10 $ exp `compose` sin (pi/2) :: [Double]
[2.718281828459045,1.664412599305428e-16,-1.3591409142295225,-1.109608399536952e-16,0.45304697140984085,4.299732548205689e-17,-0.11703713428087555,-1.2516118554300737e-17,2.5551309845882386e-2,3.0070240853853573e-18]

Computes the (truncated) Taylor series for (exp . sin) at pi/2.
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe


Re: [Haskell-cafe] Hit a wall with the type system

2007-11-29 Thread Dan Weston
I must be missing something, because to me the contract seems to be much 
simpler to express (than the Functor + Isomorphism route you seem to me 
to be heading towards):


diff :: (Eq x,
 Dense x,
 Subtractible x,
 Subtractible y,
 Divisible y x yOverX) = (x - y) - x - yOverX

class Densea where addEpsilon :: a - a
class Subtractible a where takeAway   :: a - a - a
class Divisiblea b c | a b - c
 where divide :: a - b - c

Then diff is just:

diff f a = if dx == dx' then error Zero denom else dydx
  where a'   = addEpsilon  a
dx   = a'   `takeAway`   a
dx'  = a`takeAway`   a'
dy   = f a' `takeAway` f a
dydx = dy   `divide`   dx

With the instances:

instance Dense Double where addEpsilon x = x * 1.001 + 1.0e-10
instance Dense Intwhere addEpsilon x = x + 1

instance Subtractible Double where takeAway x y = x - y
instance Subtractible Intwhere takeAway x y = x - y

instance Divisible Double Double Double where divide x y = x / y
instance Divisible IntIntIntwhere divide x y = x `div` y

and throwing in {-# OPTIONS_GHC -fglasgow-exts #-}, I get:

*Diff diff sin (0.0::Double)
1.0
*Diff diff (\x - x*7+5) (4::Int)
7

Dan Weston

Chris Smith wrote:

Hi Dan, thanks for answering.

Dan Piponi wrote:

When you specify that a function has type a - b, you're entering into
a bargain. You're saying that whatever object of type a you pass in,
you'll get a type b object back. a - b is a statement of that
contract.


Sure, that much makes perfect sense, and we agree on it.


Now your automatic differentiation code can't differentiate
any old function. It can only differentiate functions built out of the
finite set of primitives you've implemented (and other polymorphic
enough functions).


Sure.  To be more specific, here's the contract I would really like.

1. You need to pass in a polymorphic function a - a, where a is, at 
*most*, restricted to being an instance of Floating.  This part I can 
already express via rank-N types.  For example, the diffFloating 
function in my original post enforces this part of the contract.


2. I can give you back the derivative, of any type b - b, so long as b 
is an instance of Num, and b can be generalized to the type a from 
condition 1.  It's that last part that I can't seem to express, without 
introducing this meaningless type called AD.


There need be no type called AD involved in this contract at all.  
Indeed, the only role that AD plays in this whole exercise is to be an 
artifact of the implementation I've chosen.  I could change my 
implementation; I could use Jerzy's implementation with a lazy infinite 
list of nth-order derivatives... or perhaps I could implement all the 
operations of Floating and its superclasses as data constructors, get 
back an expression tree, and launch Mathematica via unsafePerformIO to 
calculate its derivative symbolically, and then return a function that 
interprets the result.  And who knows what other implementations I can 
come up with?  In other words, the type AD is not actually related to 
the task at hand.


[...]


Luckily, there's a nice way to
express this. We can just say diff :: (AD a - AD a) - a - a. So AD
needs to be exported. It's an essential part of the language for
expressing your bargain, and I think it *is* the Right Answer, and an
elegant and compact way to express a difficult contract.


Really?  If you really mean you think it's the right answer, then I'll 
have to back up and try to understand how.  It seems pretty clear to me 
that it breaks abstraction in a way that is really rather dangerous.


If you mean you think it's good enough, then yes, I pretty much have 
conluded it's at least the best that's going to happen; I'm just not 
entirely comfortable with it.





___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe


[Haskell-cafe] Hit a wall with the type system

2007-11-28 Thread Chris Smith
I was talking to a few people about this on #haskell, and it was 
suggested I ask here.  I should say that I'm playing around here; don't 
mistake this for an urgent request or a serious problem.

Suppose I wanted to implement automatic differentiation of simple 
functions on real numbers; then I'd take the operations from Num, 
Fractional, and Floating, and define how to perform them on pairs of 
values and their differentials, and then I'd write a differentiate 
function... but finding an appropriate type for that function seems to 
be a challenge.

We have:

1. Differentiating a function of the most general type (Num a = a - a) 
should produce a result of type (Num a = a - a).

2. Differentiating a function of the more specific type (Fractional a = 
a - a) should produce a result of that type (Fractional a = a - a).

3. Differentiating a function of the most specific type (Floating a = a 
- a) should produce a result of type (Floating a = a - a).

4. BUT, differentiating a function which is of a more specific type than 
(Floating a = a - a) is not, in general, possible.

So differentiate should have type A a = (forall b. A b = b - b) - a 
- a, but ONLY if the type class A is a superclass of Floating.

Two partial solutions are: I can just define the differentiate function 
for Floating; but that means if I differentiate (\x - x + 1), the 
result is a function only on floating point numbers, which is less than 
desirable.  Or, I can define several functions: say, diffNum, 
diffFractional, and diffFloating... all of which have precisely the same 
implementation, but different types and require copy/paste to make them 
work.

Any thoughts?

For reference, here's the code I kludged together.  (Again, I'm only 
playing around... so I wrote this very quickly and may have gotten some 
things wrong; don't use my code without checking it first!  In 
particular, I know that this code produces derivative functions whose 
domain is too large.)

 data AD a = AD a a deriving Eq
 
 instance Show a = Show (AD a) where
 show (AD x e) = show x ++  +  ++ show e ++  eps
 
 instance Num a = Num (AD a) where
 (AD x e) + (AD y f)   = AD (x + y) (e + f)
 (AD x e) - (AD y f)   = AD (x - y) (e - f)
 (AD x e) * (AD y f)   = AD (x * y) (e * y + x * f)
 negate (AD x e)   = AD (negate x)  (negate e)
 abs (AD 0 _)  = error not differentiable: |0|
 abs (AD x e)  = AD (abs x) (e * signum x)
 signum (AD 0 e)   = error not differentiable: signum(0)
 signum (AD x e)   = AD (signum x)  0
 fromInteger i = AD (fromInteger i) 0
 
 instance Fractional a = Fractional (AD a) where
 (AD x e) / (AD y f)   = AD (x / y) ((e * y - x * f) / (y * y))
 recip (AD x e)= AD (1 / x) ((-e) / (x * x))
 fromRational x= AD (fromRational x) 0
 
 instance Floating a = Floating (AD a) where
 pi= AD pi0
 exp (AD x e)  = AD (exp x)   (e * exp x)
 sqrt (AD x e) = AD (sqrt x)  (e / (2 * sqrt x))
 log (AD x e)  = AD (log x)   (e / x)
 (AD x e) ** (AD y f)  = AD (x ** y)  (e * y * (x ** (y-1)) +
   f * (x ** y) * log x)
 sin (AD x e)  = AD (sin x)   (e * cos x)
 cos (AD x e)  = AD (cos x)   (-e * sin x)
 asin (AD x e) = AD (asin x)  (e / sqrt (1 - x ** 2))
 acos (AD x e) = AD (acos x)  (-e / sqrt (1 - x ** 2))
 atan (AD x e) = AD (atan x)  (e / (1 + x ** 2))
 sinh (AD x e) = AD (sinh x)  (e * cosh x)
 cosh (AD x e) = AD (cosh x)  (e * sinh x)
 asinh (AD x e)= AD (asinh x) (e / sqrt (x^2 + 1))
 acosh (AD x e)= AD (acosh x) (e / sqrt (x^2 - 1))
 atanh (AD x e)= AD (atanh x) (e / (1 - x^2))
 
 diffNum:: Num b= (forall a. Num a= a - a) - b - b
 diffFractional :: Fractional b = (forall a. Fractional a = a - a) - b - b
 diffFloating   :: Floating b   = (forall a. Floating a   = a - a) - b - b
 
 diffNum f x= let AD y dy = f (AD x 1) in dy
 diffFractional f x = let AD y dy = f (AD x 1) in dy
 diffFloating f x   = let AD y dy = f (AD x 1) in dy

-- 
Chris Smith

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe


Re: [Haskell-cafe] Hit a wall with the type system

2007-11-28 Thread Luke Palmer
On Nov 29, 2007 4:02 AM, Chris Smith [EMAIL PROTECTED] wrote:
 I was talking to a few people about this on #haskell, and it was
 suggested I ask here.  I should say that I'm playing around here; don't
 mistake this for an urgent request or a serious problem.

 Suppose I wanted to implement automatic differentiation of simple
 functions on real numbers; then I'd take the operations from Num,
 Fractional, and Floating, and define how to perform them on pairs of
 values and their differentials, and then I'd write a differentiate
 function... but finding an appropriate type for that function seems to
 be a challenge.

 We have:

 1. Differentiating a function of the most general type (Num a = a - a)
 should produce a result of type (Num a = a - a).

I don't see why this should be true.  Int - Int is an instance of this type,
but derivatives require limits, which integers don't have.  Do you intend to
output the difference sequence of the function in this case?

But then Double - Double is also an instance of this type.  Do you intend
to approximate the real derivative when it's specialized to this?

Instead of worrying about the types, first tell us what semantics you want.

Luke
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe


Re: [Haskell-cafe] Hit a wall with the type system

2007-11-28 Thread Luke Palmer
On Nov 29, 2007 4:31 AM, Luke Palmer [EMAIL PROTECTED] wrote:
 On Nov 29, 2007 4:02 AM, Chris Smith [EMAIL PROTECTED] wrote:
  I was talking to a few people about this on #haskell, and it was
  suggested I ask here.  I should say that I'm playing around here; don't
  mistake this for an urgent request or a serious problem.
 
  Suppose I wanted to implement automatic differentiation of simple
  functions on real numbers; then I'd take the operations from Num,
  Fractional, and Floating, and define how to perform them on pairs of
  values and their differentials, and then I'd write a differentiate
  function... but finding an appropriate type for that function seems to
  be a challenge.

Oh, I think I totally missed the point.  I missed the word simple.

I think the problem is that a function of type Num a = a - a can be any
function whatsoever, it does not have to be a simple combination of operators
(it could, for example, use show, do a string transformation, and then read
the result). So while you can do your AD type, I think a function which
differentiates (Num a = a - a) is not possible using this approach.  You
must resort to numerical methods...

Luke

  We have:
 
  1. Differentiating a function of the most general type (Num a = a - a)
  should produce a result of type (Num a = a - a).

 I don't see why this should be true.  Int - Int is an instance of this type,
 but derivatives require limits, which integers don't have.  Do you intend to
 output the difference sequence of the function in this case?

 But then Double - Double is also an instance of this type.  Do you intend
 to approximate the real derivative when it's specialized to this?

 Instead of worrying about the types, first tell us what semantics you want.

 Luke

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe