CODEWITHSUNDEEP
haskellfunctional-programmingnumerical-methodsautomatic-differentiation
Hermann
Hermann

Reputation: 55

How to implement Carleman Matrix in Haskell?

I am trying to implement Carlemann matrix of a differentiable function in Haskell using the Numeric.AD library. I'm using https://en.wikipedia.org/wiki/Carleman_matrix for reference.

So far I have the following,

{-# LANGUAGE ParallelListComp #-}

module Sqrt where
import Data.List
import Numeric.AD

-- Find the kth order derivative
kdiff :: (a -> a) -> Int -> a
kdiff f 0 = diff f 0
kdiff f k = kdiff (diff f) (k-1)

-- Takes a function f and generates the Carleman matrix of f
carleman :: Floating a => (a -> a) -> Int ->  [[Float]]
carleman f size = [[ aux j k | k <- [0..size-1]] | j <- [0..size-1]]
    where
        aux :: Int -> Int -> Float
        aux j k = fromIntegral(kdiff (\x -> f x ^ j) k) / fromIntegral(product [1..k])

    • Couldn't match type ‘a’
                     with ‘AD s (Numeric.AD.Internal.Forward.Forward a)’
      Expected: AD s (Numeric.AD.Internal.Forward.Forward a)
                -> AD s (Numeric.AD.Internal.Forward.Forward a)
        Actual: a -> a
      ‘a’ is a rigid type variable bound by
        the type signature for:
          kdiff :: forall a. (a -> a) -> Int -> a
        at Sqrt.hs:8:1-29
    • In the first argument of ‘diff’, namely ‘f’
      In the first argument of ‘kdiff’, namely ‘(diff f)’
      In the expression: kdiff (diff f) (k - 1)
    • Relevant bindings include
        f :: a -> a (bound at Sqrt.hs:10:7)
        kdiff :: (a -> a) -> Int -> a (bound at Sqrt.hs:9:1)
   |
10 | kdiff f k = kdiff (diff f) (k-1)
   |                         ^
Failed, no modules loaded.

Edit: I've found that having a separate function for f helps.

import Numeric.AD

f :: Fractional a => a -> a
f x = x^2

-- Takes a function f and generates the Carleman matrix of f
carleman :: Int ->  [[Double]]
carleman size = [[ aux j k | k <- [0..size-1]] | j <- [0..size-1]]
    where
    aux j k = if k < length (diffs (\x -> f x ^ j) 0) 
              then (diffs (\x -> f x ^ j) 0 !! k) / product [1.. fromIntegral k]
              else 0

Upvotes: 2

Views: 84

Answers (1)

leftaroundabout
leftaroundabout

Reputation: 120751

The signature

carleman :: Floating a => (a -> a) -> Int -> [[Float]]

...which is really shorthand notation for

{-# LANGUAGE ExplicitForAll, UnicodeSyntax #-}

carleman :: ∀ a . Floating a => (a -> a) -> Int -> [[Float]]

expresses that carleman should be able to work with functions of whatever type. For example, I should be able to use it with

f :: Double -> Double
f x = fromIntegral (denominator (toRational x))

which is massively discontinuous, never mind differentiable.

To be able to differentiate the function, it can't be simply Double or Float. In this case it needs to operate, as the error message says, on AD s (Numeric.AD.Internal.Forward.Forward a) values, which take care of the additional bookkeeping involved in automatic differentiation. Actually that's not sufficient because you need higher derivatives; multiple differentiation is supported by AD s (Tower a). And the s type variable needs to be universally quantified within the function to be differentated, which requires {-# LANGUAGE RankNTypes #-}. (This essentially ensures you're not mixing up the infinitesimals of different functions that are being differentiated.)

All the tedious computing of each k-th derivative is not only a problem because each of them would need to have a different type (you need different evaluation variables for each order), and massively inefficient because you're differentiating each order separately instead of re-using the work, but also completely unnecessary because you can simply obtain the whole list of derivatives in one go and use that:

carleman :: ∀ a . Floating a
                    => (∀ s . AD s (Tower a) -> AD s (Tower a))
                    -> Int -> [[a]]
carleman f size = [ [ f' / product [1..k]
                    | (k,f') <- zip [0..size-1] (diffs (\x -> f x ^ j) 0) ]
                  | j <- [0..size-1]
                  ]

Furthermore I would point out that there's no good reason to restrict the size upfront. You can simply create an infinite matrix and trim that later on to whatever finite size you may need. And computing product [1..k] is a bad way of getting at the factorials; again it's simpler and more efficient to compute the list of all of them in one pass and go through that, then it's just a zip:

carleman :: ∀ a . Floating a => (∀ s . AD s (Tower a) -> AD s (Tower a)) -> [[a]]
carleman f = [ zipWith (/) (diffs (\x -> f x ^ j) 0) factorials
             | j <- [0..]
             ]
 where factorials = scanl (*) 1 [1..]

(I'd still be concerned about the large numbers involved; usually it's possible to reformulate such formulas in a way that avoids having large powers and factorials, but this is off topic here.)

If you find it awkward to have such an obscure type as AD s (Tower a) in your signature, an alternative is to just make the input a general numerical function, i.e. a polymorphic Floating function, and have it internally specialise to the particular differentiable type:

carleman :: ∀ a . Floating a => (∀ b . Floating b => b -> b) -> [[a]]

With either of those signatures it is possible to call the function with something like algebraic lambda expressions, e.g. carleman (\t -> sin t + t^2).

What's simply not possible though is to have the function have the same type as the resulting matrix entries. It can't be possible because extra information needs to be passed through the function in order to compute those entries.

Upvotes: 1

Related Questions